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Requirements  for  the  Degree  of  Doctor  of  Philosophy 

A DISTANCE  TO  THE  HYADES  USING  A LEAST 
SQUARES  WITH  PROBABILISTIC  CONSTRAINTS 

By 

WILLIAM  J.  COOKE 
December,  1993 

Chairman:  Heinrich  Eichhom 
Major  Department:  Astronomy 

Due  to  its  proximity,  the  Hyades  open  cluster  plays  a vital  role  in  the  establishment 
of  the  galactic  and  extragalactic  distance  scales,  as  many  of  the  distance  indicators  used  to 
estimate  the  distances  to  other  galaxies  are  calibrated  using  stars  which  are  members  of  the 
Hyades.  This  importance  of  the  Hyades  stars  as  calibrators  demands  that  their  distances 
be  determined  as  accurately  as  possible.  Past  investigations  have  made  use  of  convergent 
point  approaches,  which  impose  the  physically-ridiculous  requirement  that  all  stars  which 
are  members  of  the  cluster  must  share  exactly  the  same  velocity.  In  addition,  these 
investigations  treated  the  proper  motion  and  radial  velocity  data  in  a piecemeal  manner, 
and  have  totally  ignored  the  information  provided  by  the  trigonometric  parallaxes. 

For  the  first  time,  the  method  of  least  squares  with  probabilistic  constraints  has 
been  applied  to  the  Hyades  cluster.  This  technique  constrains  the  stars'  velocities  to  obey 
a univariate  Gaussian  (i.e.,  Maxwellian)  distribution  with  a dispersion  of  -0.25  km  s1, 
which  seems  to  closely  approximate  the  actual  velocity  distribution  of  the  Hyades,  and 
also  restricts  the  stellar  parallaxes  to  prevent  a distortion  of  the  cluster  along  the  line  of 
sight.  All  kinematic  observations  — proper  motions,  radial  velocities  and  trigonometric 
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parallaxes  — can  be  simultaneously  incorporated  into  a single  adjustment  which  directly 
yields  accurate  estimates  of  the  parallaxes  and  velocities  of  the  Hyades  stars. 

Application  of  the  technique  to  model  clusters  shows  that  the  least  squares  with 
probabilistic  constraints  does  indeed  yield  more  accurate  distances  than  those  obtained  by 
the  classical  convergent  point.  The  models  also  indicate  that  the  parameter  estimates  are 
not  terribly  sensitive  to  the  form  of  the  cluster  velocity  distribution,  the  results  being 
similar  for  both  Gaussian  and  uniform  spatial  and  velocity  distributions.  Finally,  an 
adjustment  incorporating  observations  for  140  known  members  of  the  Hyades  yields  a 
cluster  distance  of  49.7  ±1.0  parsecs,  in  good  agreement  with  recent  determinations  by 
other  investigators. 
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CHAPTER  1 
INTRODUCTION 


The  complexity  of  astrometric  instrumentation  has  increased  greatly  since  the  turn  of 
the  century.  Visual  observations  with  long  focus  telescopes  and  transit  circles  have  been 
largely  replaced  by  photographic  measurements  obtained  by  automated  measuring  engines, 
and  more  recently,  by  ground  and  space-based  observations  made  with  photoelectric 
devices  and  CCDs.  These  improved  instruments  have  made  possible  milli-arcsecond 
precision  (and,  we  hope,  a corresponding  improvement  in  accuracy)  in  proper  motion  and 
parallax  observations,  an  order  of  magnitude  better  than  the  best  observations  available  in 
the  early  1900's.  An  even  greater  improvement  has  been  made  in  the  precision  of  radial 
velocity  determinations.  Spectrometers  placed  at  the  coude  focus  of  large  reflectors  can 
now  produce  radial  velocities  with  a precision  corresponding  to  a standard  deviation  of 
0.25  km  s'1  (Griffin  et  al.,  1988),  a decrease  in  standard  deviation  by  a factor  of  more  than 
20  as  compared  to  spectrographic  observations. 

A similar  sort  of  advancement  has  occurred  in  the  mathematics  of  data  adjustment, 
which  has  its  origins  in  the  techniques  developed  by  19th  century  astronomers  to  cope 
with  the  problems  involved  in  orbit  determination.  Mathematicians  have  spent  great 
amounts  of  time  laboring  to  generalize  and  extend  these  techniques,  especially  with  regard 
to  the  method  of  least  squares,  the  adjustment  procedure  most  commonly  used. 
Unfortunately,  these  more  recent  developments  have  been  largely  ignored  by  astronomers, 
the  vast  majority  of  whom  still  employ  (often  erroneously)  the  rather  restrictive  "classical" 


1 


2 


least  squares  in  their  data  adjustments1 . The  motivation  behind  the  construction  of  more 
precise  instrumentation  and  the  development  of  general,  flexible  adjustment  procedures  is 
the  same,  namely,  to  obtain  accurate  estimates  of  certain  parameters,  and  the  use  of  both 
is  necessary  to  obtain  the  most  accurate  possible  parameter  estimates.  Today's 
astronomers  spend  huge  amounts  of  time  building  and  utilizing  marvelously  sophisticated 
instruments,  only  to  then  reduce  the  data  by  numerical  techniques  developed  long  before 
the  introduction  of  computers. 

Several  circumstances  may  help  to  explain  this  state  of  affairs.  The  first  is  that  the 
typical  astronomer  is  neither  a mathematician  nor  a statistician;  most  have  little  interest  in 
expending  the  time  required  to  learn  new  adjustment  techniques  beyond  those  learned 
while  they  were  still  students,  especially  if  those  techniques  are  complex.  Such  interest, 
where  originally  present,  is  sometimes  extinguished  by  the  language  and  notation  of  the 
mathematicians,  which  can  be  obscure  and  confusing  to  a scientist,  especially  if  there  are 
no  examples  of  the  technique's  implementation.  Another  reason,  perhaps  the  most 
prevalent,  is  that  more  general  adjustment  techniques  often  yield  parameters  little  different 
from  those  obtained  via  classical  least  squares.  Lack  of  fluency  in  programming 
languages,  unfamiliarity  with  sophisticated  modem  computers  and  the  considerable  time 
required  to  produce  complex  programs  are  obstacles  faced  by  many  astronomers  in 
implementing  modem  adjustment  techniques.  It  is  usually  far  easier  to  make  use  of  one  of 
the  many  readily  available  classical  least-squares  adjustment  routines,  thereby  reducing 
programming  time  and  the  number  of  headaches.  Eichhom  summarizes  the  situation  quite 
nicely: 

Unfortunately,  most  users  of  the  least-squares  adjustment  technique  regard  it  as  a 
mere  trade  tool,  which  comes  pre-packaged  and  in  various  degrees  of  sophistication. 
If  the  tool  is  complicated,  learning  how  to  use  it  may  consume  a lot  of  time  and  in 
addition  be  difficult,  especially  when  — as  happens  not  infrequently  — the  author  of 


1 In  "classical"  least  squares  only  one  of  the  observations  (which  are  assumed  to  be  uncorrelated)  appears 
in  each  condition  equation,  which  must  be  linear. 
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the  instructions  for  using  the  tool  is  less  than  perfectly  lucid.  One  also  encounters 
often  the  — much  wiser  than  immediately  apparent  — attitude  that  "if  it  ain't  broken, 
don't  fix  it!",  which  produces  a reluctance  in  much  of  the  community  of  users  to 
abandon  a tried-and-true,  even  though  not  optimally  accurate  adjustment 
technique. . . (Eichhom,  1993,  p.  338) 

Such  is  the  case  with  past  determinations  of  the  distances  to  the  stars  which  are 
considered  to  be  physical  members  of  the  Hyades  star  cluster.  This  cluster's  close 
proximity  and  importance  as  a calibrator  in  the  galactic  and  extragalactic  distance  scales 
have  stimulated  many  kinematic  investigations,  beginning  with  that  of  Boss  (1908).  The 
precision  of  the  observations  has  improved,  but  alas,  the  sophistication  of  the  adjustment 
techniques  used  to  determine  the  distances  to  the  member  stars  has  not.  All  past 
determinations  have  neglected  to  incorporate  the  various  types  of  available  kinematic  data 
(proper  motions,  radial  velocities  and  parallaxes)  simultaneously  into  a single  adjustment, 
opting  instead  to  treat  the  data  in  a piecemeal  fashion.  Indeed,  the  observed  parallaxes,  if 
considered  at  all,  are  simply  averaged  using  some  sort  of  weighting  scheme  and  used  as  a 
"check"  on  the  results  of  the  adjustment  procedure  in  which  the  parallax  has  been 
eliminated  from  the  condition  equations.  Further,  the  nonlinear  equations  of  condition  are 
often  linearized  in  such  a manner  as  to  introduce  systematic  errors  (i.e.,  biases)  in  the 
estimated  parameters,  especially  in  the  investigations  performed  in  the  first  decades  of  this 
century.  There  was  no  realistic  alternative  to  using  such  methodology  in  investigations 
conducted  before  high  speed  computers  became  readily  accessible,  as  an  adjustment 
involving  observations  of  a few  hundred  stars  and  multiple  equations  of  condition  can 
strain  the  capabilities  of  even  today's  computers.  Yet,  the  most  recent  published 
investigations  still  fail  even  to  incorporate  the  parallax  observations  into  the  adjustments. 

Previous  kinematic  distance  determinations  have  assumed  that  all  stars  belonging  to 
the  Hyades  cluster  move  with  exactly  the  same  velocity,  which  makes  no  physical  sense. 
As  Hanson  (1975)  points  out,  the  "principle  would  be  inapplicable  in  the  case  of  the 
Hyades  if . . . random  motions  of  individual  members  contributed  significantly  to  the  stars' 
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space  velocities"  (p.  393).  The  assumption  of  a common  space  motion  has  yielded 
acceptable  results  only  because  open  clusters  like  the  Hyades  are  loosely  bound,  with 
velocity  components  dispersed  on  the  order  of  0.5  km  s*1.  This  is  about  1%  of  the  Hyades 
space  speed,  which  is  approximately  44  km  s'*.  Yet,  adjustment  techniques  exist  that 
would  allow  the  individual  stellar  velocities  to  vary  according  to  a Gaussian  distribution 
(Eichhom,  1977).  While  it  is  probable  that  the  velocity  distribution  function  of  the  Hyades 
stars  may  not  be  Gaussian  in  nature,  an  approach  of  this  sort  is  more  physically  realistic 
than  imposing  the  physically  impossible  and  rigorous  constraint  that  all  stars  share  the 
same  velocity. 

These  same  investigations  determined  an  individual  stellar  parallax  by  placing  the  star 
at  a distance  such  that  the  space  velocity  of  the  star  corresponded  to  its  observed  proper 
motion.  However,  the  unavoidable  errors  in  the  proper  motions  will  thus  cause  the  cluster 
to  appear  more  elongated  along  the  line-of-sight  direction  than  it  actually  is.  Again,  more 
accurate  estimates  of  the  stellar  parallaxes  may  be  obtained  by  assuming  a spherical  shape 
for  the  cluster  and  stochastically  constraining  the  stars  to  lie  within  these  bounds 
(Eichhom,  1977). 

The  above  paragraphs  summarize  the  stimuli  for  the  work  presented  in  this  paper, 
which  has  two  primary  purposes:  The  first  of  these  is  the  goal  of  any  adjustment,  namely, 
to  obtain  from  the  observed  data  the  most  accurate  available  estimates  of  the  unknown 
parameters,  which  in  this  case  are  the  velocities  and  parallaxes  of  the  Hyades  stars.  The 
second,  and  more  important  aim,  is  to  break  with  the  past  and  apply  for  the  first  time  a 
modem,  general  least-squares  adjustment  technique  to  the  Hyades  kinematic  observations. 
It  is  hoped  that  this  investigation  will  be  a step  toward  altering  the  perceptions  of  the 
astronomical  community  regarding  general  least-squares  adjustments.  Like  modem 
instruments,  these  techniques  are  complex  and  sometimes  problematic,  but  only  they  can 
fully  utilize  the  existing  data  to  produce  the  most  accurate  results.  Today's  observations 


deserve  no  lesser  treatment. 


CHAPTER  2 

PREVIOUS  HYADES  DISTANCE  DETERMINATIONS 

Overview 

Past  distance  determinations  to  the  Hyades  stars  may  be  divided  into  two  broad 
groups:  Those  making  use  of  geometric-kinematic  data  (proper  motions,  radial  velocities 
and  parallaxes)  and  those  obtained  by  astrophysical  considerations  (photometric 
"parallaxes",  stellar  models,  the  Wilson-Bappu  effect,  etc.).  An  examination  of  the 
literature  also  suggests  that  the  date  of  an  investigation  generally  places  it  into  one  of 
three  distinct  "epochs,"  each  characterized  by  a particular  school  of  thought  and  use  of 
adjustment  procedures.  These  epochs  are  summarized  in  the  paragraphs  below. 

Historic  (1908-1945).  This  epoch  commences  with  the  pioneering  convergent  point 
solution  by  Boss  (1908),  followed  in  short  order  by  investigations  conducted  by  such 
notables  as  Plummer  (1912),  Rasmuson  (1921),  and  Smart  (1939).  These  investigations 
are  characterized  by  convergent  point  solutions  using  Charlier's  method  (Charlier,  1916) 
or  a variant  thereof,  and  utilize  proper  motions  on  the  system  of  the  General  Catalogue 
(GC)  with  sparse  radial  velocities  of  (by  contemporary  standards)  poor  precision.  The 
distance  to  the  Hyades1  was  determined  to  be  under  40  parsecs,  a representative  value 
being  that  of  36.3  ± 1.2  parsecs2  obtained  by  Smart.  This  period  was  brought  to  an 
abrupt  end  by  Frederick  Seares  (1945),  who  pointed  out  that  the  coefficients  of  the 
unknown  parameters  in  the  Charlier  equations  of  condition  involved  the  proper  motion 
observations  and  that  the  errors  in  these  observations  would  introduce  a bias  into  the 

1 By  "the  distance  to  the  Hyades",  the  author  refers  to  the  distance  of  the  geometrical  center  of  the  cluster 
as  determined  by  the  stars  considered  in  the  adjustment. 

2 The  numbers  given  after  the  ± are  standard  errors,  unless  otherwise  indicated. 
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convergent  point  solution.  He  developed  a new  set  of  normal  equations  which  were 
applied  to  Smart's  material,  resulting  in  a cluster  distance  of  39.2  ±1.2  parsecs,  a change 
of  some  7%. 

Classic  (1950-1966).  This  epoch  may  also  be  referred  to  as  the  "heyday  of  the 
convergent  point,"  for  it  is  during  this  period  that  the  "definitive"  investigations  of  Brown 
(1950),  Van  Bueren  (1952)  and  Wayman  et  al.  (1965)  were  conducted.  They  utilized 
transit  circle  proper  motions  either  on  the  system  defined  by  the  GC  (Van  Bueren,  1952) 
or  that  of  the  N30  (Wayman  et  al.,  1965)  and  more  numerous  and  precise  radial  velocities 
from  Wilson's  catalog.  In  order  to  avoid  the  bias  associated  with  early  convergent  point 
solutions,  a new  technique  was  introduced  in  which  an  approximation  to  the  position  of 
the  convergent  point  was  assumed  (usually  that  determined  by  Smart  (1939)),  which  was 
then  corrected  by  a least-squares  adjustment  performed  only  with  the  use  of  the  position 
angles  of  the  observed  proper  motions.  Though  this  technique  proved  to  be  successful  in 
avoiding  any  obvious  bias,  information  was  lost  as  the  magnitudes  of  the  proper  motions 
were  removed  from  the  adjustment  procedure.  Wayman  et  al.  obtained  a Hyades  distance 
of  40.6  ±1.6  parsecs,  in  good  agreement  with  the  value  of  40.4  ±1.1  parsecs  determined 
by  Van  Bueren.  The  good  agreement  between  the  convergent  point  solutions  led  many 
astronomers  to  feel  that  the  question  of  the  Hyades  distance  was  closed  and  that  more 
work  would  yield  only  very  minor  corrections,  hardly  worth  the  effort  involved.  This 
feeling  was  shattered  by  the  controversial  paper  of  Hodge  and  Wallerstein  (1966),  who 
argued  that  the  Hyades  distance  ought  to  be  increased  to  48.4  parsecs  based  on  (a)  the 
mean  of  the  parallaxes  of  Hyades  stars  in  the  Yale  parallax  catalog,  (b)  the  dynamical 
parallaxes  of  Hyades  visual  binaries  and  (c)  the  Wilson-Bappu  effect. 

Modem  (1967-present).  The  current  epoch  is  marked  by  high  precision  photographic 
proper  motions  and  parallaxes  obtained  through  the  use  of  automated  measuring 
machines,  precise  radial  velocities  obtained  by  spectrometers,  and  the  widespread 
availability  of  computers.  The  convergent  point  solutions  have  been  recast  into  new 
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forms:  the  proper  motion  gradient  method  employed  by  Upton  (1970),  and  the  radial 
velocity  gradient  method  recently  developed  by  Griffin  et  al.  (1988).  Even  though  the 
direct  determination  of  a convergent  point  is  avoided  by  these  techniques,  the  assumption 
of  a common  cluster  motion  is  maintained  and  the  stellar  parallaxes  are  still  obtained  in  a 
manner  similar  to  the  old  convergent  point  solutions.  In  addition,  a new  emphasis  has 
developed  on  cluster  distance  determined  via  "astrophysical"  or  "secondary"  means,  such 
as  R-I  photometry.  These  techniques  produce  results  which  are  accorded  a good  deal  of 
weight  and  are  used  as  comparisons  for  the  convergent  point  solutions.  Trigonometric 
parallaxes  continue  to  be  left  out  of  the  still-primitive  adjustment  procedures,  the  mean 
sometimes  being  calculated  for  comparison  with  the  results  of  an  adjustment.  Typical 
modem  values  for  the  Hyades'  distance  are  47.2  ± 1.9  parsecs  (convergent  point  solution, 
Hanson,  1980),  45.3  ±1.6  parsecs  (mean  of  trigonometric  parallaxes,  Upgren,  1974a)  and 
44.0  ± 0.4  parsecs  (R-I  photometry,  Upgren,  1974b). 

Convergent  Point  Determinations 


Derivation  of  the  Classical  Convergent  Point  Condition  Equations 


Let  x be  the  location  of  a star  which  is  a member  of  the  cluster  and  x be  its  velocity. 
Then  we  may  write 


^cosacosb^ 


x = rx(a,90°-8)  = r 


sin  a cos  5 
sin8 


(2-1) 


x = rx  + rx 


where  a,  8 are  the  star's  right  ascension  and  declination,  r is  its  distance,  and  r its  the 
radial  velocity,  x , the  time-derivative  of  the  star's  direction  unit  vector,  is  given  by 


r x 


8 


O 


Figure  2-1.  Convergent  point  geometry 


X = 


^-dcos8sina-8sin8cosa^ 

dcos8cosa-8sin8sina 

8cos8 


(2-2) 


Note  that  the  units  of  & and  8 are  radians  per  unit  time,  whereas  the  commonly 
communicated  |ia  and  |i8  are  in  units  of " yr1  or  " per  century.  A little  algebraic 

manipulation  enables  one  to  show  that  the  magnitude  of  x is  the  total  proper  motion  of 
the  star,  (i: 


Ixl  = Vex2  cos28+82  = it 


(2-3) 


By  definition,  x • x = 0 . It  can  also  be  shown  that 

r = xx, 

so  that  the  second  of  equations  (2-1)  becomes 

x = (x  • x)x  + rx , or  rx  = x - (x  • x)x 


Now  form 


rx  ■ rx  = (x  - (x  • x)x)  • (x  - (x  • x)x), 


(2-4) 


(2-5) 


(2-6) 


which  is  equivalent  to 

r2|x|2  = |x| 2 -r2,  or  r2p.2  =V2  -r2  (2-7) 

where  V is  the  magnitude  of  the  star's  space  velocity  (expressed  in  units  of  km  s1),  i.e., 
the  star's  speed.  A glance  at  figure  2-1  gives  an  alternative  expression  for  the  radial 
velocity: 


r = V cos  A,, 


(2-8) 


10 


where  A,  is  the  angle  between  the  line-of-sight  direction  and  the  star's  velocity. 

Substituting  the  above  into  equation  (2-7)  results  in 

r2\i2  -V2  - V2  cos2  A,  = F2(l-cos2  a)  = V2  sin2  A.,  or 

V sin  A,  (2-9) 

r = 

kh 


where  K is  the  appropriate  factor  for  converting  the  radial  velocities  from  units  of  A.U.  yr1 
to  units  of  km  s1  (~  4.74).  In  terms  of  the  parallax, 


Kjl 

Esin  A 


(2-10) 


Equation  (2-10)  is  a general  result,  applicable  to  any  star,  whether  or  not  it  is  a cluster 
member,  provided  that  A is  the  angle  between  the  star's  position  and  its  velocity. 

However,  in  the  case  of  a group  of  stars  sharing  a common  velocity,  the  stars  will 
appear  to  converge  to  (or  diverge  from)  a point3  in  the  sky,  called  the  convergent 
(explosion)  point.  The  direction  to  this  point  is  the  same  as  the  group's  common  space 
velocity,  and  may  be  written 


^cos  AcosD^ 
x = Vx(  A,90°-D)  = V sin  AcosD 

sinD 


V 


y 


(2-11) 


where  A and  D are  the  right  ascension  and  declination,  respectively,  of  the  convergent 
point.  Once  A and  D are  known,  A is  computed  from 


cosA  = sin5sinZ)  + cos5cosZ)cos(a-A) 


(2-12) 


3 In  actuality,  this  will  be  an  area  rather  than  a "point",  due  to  the  unavoidable  errors  in  the  measured 
proper  motions  and  the  fact  that  the  stars  do  not,  in  fact,  have  exactly  the  same  velocity. 
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A knowledge  of  each  star's  X,  combined  with  the  cluster  velocity  and  the  observed  proper 
motions,  enable  one  to  carry  out  the  computation  of  the  individual  stellar  parallaxes  by 
means  of  equation  (2-10).  Note  that  this  equation  involves  only  the  total  proper  motion 
and  not  its  direction,  and  that  it  is  independent  of  the  radial  velocity.  It  does,  however, 
require  knowledge  of  the  cluster  velocity.  As  all  stars  share  this,  the  velocity  may  be 
estimated  through  a standard  least-squares  adjustment,  with  the  velocity  components  as 
the  global  variables. 

The  equations  of  condition  are  set  up  by  inverting  the  second  of  equations  (2-1), 
which  may  be  rewritten  in  the  form  (Eichhom,  1982): 


( 


x = rR3(-a)R2(8)  acosS 

8 


7t  r 


(2-13) 


where  R2  and  R3  are  standard  rotation  matrices,  defined  by 


^cos0 

0 

-sin0" 

^ COS(J) 

sin<{) 

0^ 

R2(9)  = 

0 

1 

0 

, R3(<t>)  = 

— sin  <|> 

COS(|) 

0 

^sin0 

0 

cos0  J 

, 0 

0 

h 

We  solve  equation  (2-13)  with  respect  to  the  vector  on  the  right-hand  side  and  obtain 


f 


\ 


nr 
a cos  8 
8 


\ 


= 7tR,(-8)R2(a)x, 


(2-14) 


which  may  be  explicitly  written 


r-xcosa  cos8  -ysina  cos8  -zsinS 
a cos8  + x7t  sin  a - yn  cosa 
v8  + x7t  cosa  sin8  + y7t  sina  sin8  - Z7t  cos8 


\ 

= 0, 


(2-15) 
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where  x,  y and  z are  the  cartesian  components  of  x (units  are  km  s1)- 

These  are  the  desired  condition  equations  necessary  for  the  determination  of  the 
cluster  velocity.  Note  that  they  are  nonlinear  if  one  regards  the  parallax,  k,  as  an 
observation,  because  of  the  terms  m,  yrc  , and  zk  . Past  investigators  have  employed 
various  techniques  for  circumventing  the  failure  of  equations  (2-15)  to  be  linear,  which 
shall  be  considered  in  some  detail  in  the  following  sections. 

Charlier's  Method 

Before  1945,  investigators  employed  a method  developed  by  Charlier  (1916)  to  obtain 
the  distances  to  the  stars  of  the  Hyades  cluster.  This  technique  ignored  the  first  of 
equations  (2-15)  altogether,  wasting  practically  all  the  information  contained  therein.  The 
parallax  was  then  eliminated  from  the  second  and  third  of  equations  (2-15),  resulting  in  the 
condition  equation 


ax  + by- cz  — 0. 


(2-16) 


where 

a = acos8cosasinS-8sina 
b = acos8sinasin8  + 8cosa 
c = acos2  8 

In  order  to  simplify  the  computations,  equation  (2-16)  was  rewritten  as 

a--c-+b  = 0 (2-17) 

y y 

and  a regression  performed  to  estimate  values  for  the  ratios  of  the  velocity  components. 
The  coordinates  of  the  convergent  point  were  then  obtained  from  equations  (2-11).  Note 
that  the  quantity  minimized  in  this  "adjustment"  has  nothing  to  do  with  the  sums  of  the 
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squares  of  any  residuals.  The  only  justification  for  employing  this  wretched  model  is  that 
it  leads  to  a result  (after  a fashion),  with  the  minimum  arithmetic  effort,  an  important 
consideration  in  those  days!  V is,  of  course,  still  missing  completely;  but  equation  (2-8), 
recast  in  the  form 


V = rsecX,  (2-18) 

became  the  condition  equation  for  the  second  of  the  adjustments,  now  with  the  aim  to 
estimate  V,  involving  the  observed  radial  velocities.  The  cluster  velocity  having  thus  been 
determined  by  these  two  separate  adjustments,  the  individual  stellar  parallaxes  were 
computed  by  means  of  equation  (2-10). 

At  this  point,  one  may  ask  why  Smart  and  other  early  investigators  went  to  the  trouble 
of  eliminating  the  parallax  from  equation  (2-15),  when  they  could  obtain  the  velocity 
components  directly  from  considering  only  the  first  of  these  equations,  i.e.,  by  working 
with  the  radial  velocities  exclusively.  The  answer  is  two-fold:  one,  the  Hyades  cluster  is 
fairly  concentrated,  with  most  stars  lying  within  a few  degrees  of  the  cluster's  center,  and 
two,  the  precision  of  the  available  radial  velocities  was  rather  poor  (standard  errors 
greater  than  5 km  s'1).  Even  with  today's  high  precision  radial  velocities  (standard  errors 
of  less  than  1 km  s'1),  these  reasons  still  make  the  first  of  equations  (2-15)  a poor  choice 
for  an  adjustment  condition  equation.  Why  a simultaneous  adjustment  of  equations  (2-15) 
and  (2-16)  was  not  performed  remains,  however,  still  puzzling.  Surely  these  investigators 
were  familiar  with  the  concept  of  weight! 

Many  investigators  used  this  technique  to  estimate  distances  to  the  Hyades  stars, 
including  Plummer  (1912),  Rasmuson  (1921),  Smart  (1939),  Losert  (1939),  and  Pearce 
(1955).  All,  with  the  exception  of  Boss  (1908),  obtained  a Hyades  distance  under  40 
parsecs,  with  Smart's  value  of  36.3  ±1.2  parsecs  being  representative.  Table  2-1 
summarizes  the  results  of  some  of  these  investigations. 
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Table  2-1.  Early  Hyades  distance  determinations. 


Investigator 

a CP 

6 CP 

Distance  (pc) 

No.  stars 

Boss  (1908) 

92°.82  ± 2°.2 

6°.92  ± 0°.44 

41.7 

41 

Losert  (1939) 

92°.75  ± 1°.6 

6°.87  ± 1°.6 

37 

10 

Smart  (1939) 

91°.83±0°.66 

8°.07  ± 0°.26 

36.3  ± 1.2 

72 

Pearce  (1955) 

91°.85  ± 0°.65 

9°.l  ±0°.25 

37  ± 1 

74 

Seares'  Criticism  of  Charlier's  Method 


So  matters  stood  until  1945,  when  Frederick  Seares  published  a critique  of  Charlier's 
method,  in  which  he  showed  that  this  technique  would  introduce  a bias,  i.e.,  systematically 
underestimate  the  distances  to  Hyades  members.  The  problem  arises  because  the 
coefficients  of  condition  equations  (2-17)  contain  the  observed  and  therefore  unavoidably 
error- affected  proper  motion  components  of  the  star.  Incorporating  these  errors  into  the 
condition  equations  yields 

a'--c'-+b'  = 0 (2-20) 

y y 


where  a',  b',  and  c'are  given  by 


a'=a+  cosdcosasinS-e^  sina 

b'  = b + E^  cos8sinasin8  + e^ cosa  (2-21) 

c'=  c + E^  cos8 

with  and  being  the  errors  in  the  proper  motion  components.  The  primed 

coefficients  are  the  ones  unwittingly  used  in  Charlier's  method;  the  unprimed  coefficients 
are  the  "true"  values.  The  ordinary  least-squares  solution  of  equations  (2-20)  is  given  by 

A = -(ata)_1AtF0  (2-22) 


where  A is  the  vector  of  parameter  estimates,  and 
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At  = 


f a\ 


Vci 


-c' 


F„  = 


\ 


(2-23) 


with  m being  the  number  of  stars.  ATA  has  form 


AtA  = 


I aK  'Z-a'c' 
1=1  1=1 
m m 

I -ay,  5>' 


V >=i 


1=1 


or,  expanding  the  first  column, 


AtA  = 


m / \ iw 

yaci+  2c3£(jaA[  - 2c£,^A.2  — 2£(ja£^A1  A2  + 6^  AI+£>ig  A2j 

i=l  i=l 

m / \ m 


V i=l 


i=l  V 


where 

A,  = cosScosasinb 
A2  = sin  a 
A3  = cos  5 

Note  that  the  elements  of  ATA  differ  from  the  "true"  elements  by  terms  involving  the 
errors  in  the  proper  motion  components.  The  terms  containing  e^2  and  £^2  account  for 

most  of  the  differences,  as  the  linear  and  mixed  terms  mostly  cancel.  These  terms  will 
introduce  systematic  errors  into  the  parameter  estimates,  and  hence  the  direction  to  the 
convergent  point. 

In  applying  the  above  to  Smart's  data,  Seares  revised  the  equatorial  coordinates  of  the 
Hyades  convergent  point  from  Smart's  original  (91°.83,  8°.07)  to  (93°.48,  7°.53).  The 
alteration  in  the  convergent  point  direction  resulted  in  a change  to  the  cluster  distance  by 
some  7%,  moving  the  center  of  the  Hyades  out  by  approximately  3 parsecs  to  39.2  ± 1.2 
parsecs.  It  should  perhaps  be  emphasized  that  Seares  kept  only  the  terms  involving  E^2 
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and  e^52  in  the  normal  equations  (apparently  because  he  trusted  the  linear  — and 
Ejxg  — and  mixed  — £|iae|a8  — terms  to  cancel),  rendering  his  correction  to  the  Hyades 
distance  no  more  than  a first  order  approximation.  Nonetheless,  his  paper  caused  a minor 
stir  in  the  astrometric  community  which  resulted  in  the  development  of  a new  convergent 
point  technique. 

The  Classical  Convergent  Point 

Brown  (1950)  also  pointed  out  the  shortcomings  in  Charlier's  method,  but  was  equally 
critical  of  Seares'  regression  correction.  His  dissatisfaction  with  the  techniques  then 
employed  caused  him  to  evaluate  different  approaches  in  solving  for  the  convergent  point, 
leading  to  his  conclusion  that  it  would  be  better  to  use  the  position  angles  of  the  individual 
stellar  proper  motions  rather  than  the  proper  motion  components.  While  keeping  the 
assumption  of  identical  stellar  velocities.  Brown  proposed  the  abandonment  of  obtaining 
the  condition  equations  from  equations  (2-15)  in  favor  of  a derivation  from  the  equation 
giving  the  position  angle  at  the  star's  position  (a,  5)  of  the  great  circle  through  (a,  8)  and 
the  coordinates  of  the  convergent  point  (A,  D).  He  also  advocated  obtaining  a solution 
via  differential  corrections  or  an  alternate  maximum  likelihood  technique  rather  than  using 
a least  squares  regression.  However,  his  application  of  the  maximum  likelihood  approach 
to  the  Ursa  Major  association  met  with  very  limited  success,  as  the  numerical 
computations  were  quite  complicated  for  the  available  computational  facilities,  i.e., 
mechanical  calculators. 

Later  investigators  (van  Bueren,  1952;  Wayman  et  al.,  1965)  adopted  Brown's  concept 
of  using  the  proper  motion  position  angles  to  determine  the  convergent  point,  but  returned 
to  least  squares  as  the  technique  used  to  obtain  the  solutions.  Starting  with  the  great 
circle  equation  used  by  Brown, 


cot0  = tanDcosScsc(.4-a)-sin8cot(.4-a), 


(2-24) 
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with  0 being  the  position  angle  at  the  star's  position  (a,  8)  of  the  great  circle  through  (a, 

8)  and  the  coordinates  of  the  convergent  point  (A,  D),  they  performed  a differentiation 
with  respect  to  A and  D.  This  resulted  in  the  condition  equation  for  the  adjustment, 

BAA -CAD  = 0O  -6  , (2-25) 


where  0q  are  the  stellar  proper  motion  position  angles  (measured  from  the  direction  to  the 
pole  of  rotation),  determined  from 


acos8 
tan0o  = — t — 


(2-26) 


and  B and  C are  given  by 


B=  cot(^40  -a)cos0sin0-sin8sin2  0 
C = sec2  D0  cos8csc(4)  - oc)sin2  0 


(2-27) 


Note  that  B and  C do  not  depend  on  any  observed  quantities  (the  star  positions  are  known 
with  great  precision  and  so  are  treated  as  constants),  thereby  avoiding  the  regression  error 
associated  with  Charlier's  technique.  The  parameters  obtained  by  the  adjustment  are  AA 
and  AD,  these  being  the  corrections  to  the  assumed  position  of  the  convergent  point  (Aq, 
D0).  The  coordinates  of  the  convergent  point  were  then  found  from  the  simple  relations 


A = A0+AA 
D = D0+ AD 


(2-28) 


Once  these  were  obtained,  the  cluster  velocity  and  the  individual  stellar  parallaxes  were 
determined  as  in  Charlier's  method.  Brown,  in  conformance  with  then  and  — alas  — even 
now  prevailing  common  practice  makes  no  attempt  to  justify  the  minimization  of  the  sums 
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of  the  squares  of  the  position  angle  residuals,  although  such  a justification  could  be 
contrived  by  invoking  the  central  limit  theorem. 

The  first,  and  perhaps  best  known,  investigation  using  the  above  technique  was  that  of 
van  Bueren  (1952),  who  placed  the  Hyades  convergent  point  at  (94°.62  ± 0°.53,  7°.48  ± 
0°.18)  and  the  distance  to  the  cluster  center  at  40.4  ±1.1  parsecs.  Since  then  there  have 
been  many  investigations  of  the  Hyades  distance  using  the  classical  convergent  point 
approach,  including  a detailed  study  by  Wayman  et  al.  (1965).  The  results  of  some  of 
these  investigations  are  summarized  in  table  2-2  and  the  positions  of  the  convergent  points 
are  displayed  in  figure  2-2.  Note  the  gradual  progression  of  the  convergent  point  towards 
higher  right  ascensions  and  smaller  declinations,  due  in  part  to  the  systematic  differences 
among  the  proper  motion  sets  used  in  the  solutions.  These  differences  will  be  discussed  in 
some  detail  in  a later  section. 


Table  2-2.  Classical  Hyades  convergent  point  solutions. 


Investigator 

a CP 

5 CP 

Distance  (pc) 

van  Bueren  (1952) 

94°.62  ± 0°.53 

7°.48±0°.18 

40.4  ± 1.1 

Heckman  and  Ltibeck  (1956) 

94°.0 

7°.6 

41.2 

Wayman  et  al.  (1965) 

94°.  13  ± 0°.49 

7°.68±0°.19 

40.6  ± 1.6* 

Corbin  et  al.  (1975) 

96°.3 

7°.27 

43.55  ± 0.74 

Hanson  (1975) 

102°.5  ± 3°.7 

5°.9  ± 0°.8 

48.4  ± 4.8 

Morris  and  Luyten  (1983) 

95°.2  ± 0°.5 

7°.8  ± 0°.2 

42.7 

Schwan  (1991) 

97°.68  ± 0°.25 

5°.98±0°.18 

47.9  ± 0.8 

Gradient  Techniques 


Upton  (1970)  developed  a technique  in  which  the  gradients  in  either  the  right 
ascension  proper  motions  or  the  declination  proper  motions  could  be  combined  with  the 
average  radial  velocity  of  the  Hyades  stars  to  obtain  the  cluster's  distance.  Although  the 
method  does  not  directly  determine  the  position  of  the  convergent  point,  it  is  a 
"convergent  point"  technique  in  the  sense  that  it  retains  the  assumption  of  a common 
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space  velocity  for  the  Hyades  stars.  Hanson  notes  that "...  the  method  of  analysis  is 
simply  the  CP  principle  treated  in  different  mathematical  form"  (Hanson,  1975,  p.  395)4  , 
pointing  out  that  Upton's  equations  can  be  derived  by  differentiating  equation  (2-9). 

The  basic  concepts  behind  Upton's  approach  can  be  grasped  by  considering  the 
Hyades  to  be  spherical  in  shape,  having  some  finite  diameter,  s,  and  center  at  some 
distance,  r.  Then  simple  trigonometry  leads  to  the  relation 

0=-,  (2-29) 

r 

where  0 represents  the  angular  diameter  of  the  cluster.  Under  the  assumptions  that  the 
cluster  is  of  small  angular  extent  and  undergoes  no  expansion  or  contraction,  the 
differentiation  of  this  with  respect  to  time  results  in 


0 = — rr,  or 


0 

0 


r 

r 


(2-30) 


Simply  put,  this  equation  states  that  to  first  order,  the  distance  to  the  Hyades  can  be 
determined  by  dividing  its  mean  radial  velocity  by  the  fractional  decrease  in  the  cluster 
angular  diameter.  This  decrease  in  the  cluster's  angular  size  as  it  recedes  can  be 
determined  from  the  convergence  of  the  member  stars'  proper  motions;  indeed,  Upton 
(1970)  shows  that  equation  (2-30)  can  be  approximated  by 


L _ 

r da  db 


(2-31) 


4 This  statement  seems  somewhat  contrived.  The  simple  convergent  point  method  is  based  on  the  first 
time  derivatives  of  the  location  components  in  polar  coordinates,  while  the  gradient  techniques  regard  the 
(much  less  precise!)  second  time  derivatives  of  these  coordinates  as  the  observed  quantities,  thus 
neutralizing  constant  biases.  Insofar  as  Hanson  implies  that  both  classes  of  methods,  applied  to  the  same 
observations,  will  yield  the  same  results,  he  is  wrong. 
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where  the  terms  on  the  right  represent  the  proper  motion  gradients  in  right  ascension  and 
declination,  respectively.  It  should  be  noted  that  the  cluster  distance  can  be  obtained  from 
either  gradient,  and  much  has  been  made  of  the  notion  that  any  differences  between  the 
distance  determined  from  the  right  ascension  proper  motion  components  and  that 
determined  from  the  declination  components  should  give  some  idea  of  the  accidental  and 
systematic  errors  in  these  proper  motions  (see  Upton,  1970  and  Hanson,  1975).  While 
this  may  be  true  of  the  accidental  errors  involved,  it  seems  not  to  be  the  case  with  the 
systematic  errors,  as  Upton's  Hyades  distance  determination  fails  to  show  the  systematic 
errors  in  the  transit  circle  proper  motions  used  in  his  adjustment.  Similarly,  Hanson's 
(1975)  application  of  the  technique  failed  to  reveal  the  magnitude  dependency  in  his 
photographic  proper  motion  estimates  (McAlister,  1977). 

A variant  of  Upton's  approach,  one  using  the  radial  velocity  gradient,  was  developed 
by  Gunn  et  al.  (1988).  This  technique  obtains  its  basic  equation  by  differentiating  equation 
(2-8)  with  respect  to  X: 


— = V sinX  = r tanX 

dk 


(2-32) 


Substitution  of  the  above  into  equation  (2-9)  gives  the  result 


r = 


rdT 

Kit 


(2-33) 


where  p.  is  the  mean  cluster  proper  motion.  The  accurate  determination  of  the  radial 
velocity  gradient  has  only  recently  been  made  possible  by  high  precision  radial  velocities 
measured  by  Detweiler  (Detweiler  et  al.,  1984)  and  Griffin  (Griffin  et  al.,  1988),  as  the 
scatter  in  the  Wilson  radial  velocities  was  much  too  large  (of  the  order  of  5 km  s_1). 
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Figure  2-2.  Hyades  convergent  point  positions  (orientation  1950.0) 
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Table  2-3  summarizes  the  Hyades  distance  determinations  obtained  using  the  gradient 
techniques  outlined  above.  There  is  the  same  large  scatter  in  the  Hyades  distance  as  is 
seen  in  the  classical  convergent  point  investigations,  illustrating  the  difficulties  involved  in 
obtaining  absolute  proper  motions  and  radial  velocities  free  from  systematic  effects. 

Table  2-3.  Gradient  technique  distance  determinations. 

Investigator  Gradient  Type  Distance  (pc) 

Upton  (1970)  Proper  Motions  41.6  ±1.0 

Hanson  (1975)  Proper  Motions  48.4  ±4.4 

Gunn  et  al.  (1988) Radial  Velocity 45.4  ±2.1 

Secondary  Distance  Determinations 
The  Hodge  and  Wallerstein  Controversy 

In  the  mid-1960's,  Hodge  and  Wallerstein  (1966)  created  a stir  in  astronomical  ranks 
by  suggesting  that  the  distances  to  the  Hyades  as  derived  from  the  convergent  point 
determinations  were  too  small  by  some  8 parsecs.  They  based  their  argument  on  the 
results  of  investigations  involving  so-called  "secondary"  indicators  — trigonometric 
parallaxes,  dynamical  parallaxes  of  Hyades  visual  binaries  and  the  widths  of  the  K- 
emission  line  cores  in  four  Hyades  yellow  giants  (the  Wilson-Bappu  effect).  The  distances 
obtained  from  these  methods  and  presented  by  Hodge  and  Wallerstein  in  their  paper  are 
summarized  in  table  2-4. 

Although  there  was  much  initial  resistance,  it  has  now  become  accepted  that  the 
"secondary"  Hyades  distance  indicators  should  be  given  weight  equal  to  convergent  point 
techniques,  especially  with  the  evidence  mounting  that  proper  motions  are  quite  sensitive 
to  various  systematic  errors.  As  Hodge  and  Wallerstein  (1966)  point  out,  a systematic 
error  of  only  0.0008  " yr 1 deg'1  in  the  declination  proper  motions  can  cause  the  Hyades 
distance  to  change  by  almost  20%.  A glance  at  the  previous  section  also  shows  that 
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convergent  point  solutions  incorporating  modem  proper  motions  placed  on  the  FK5 
system  or  tied  to  a reference  frame  of  external  galaxies  (which  are  presumably  less 
affected  by  systematic  errors)  yield  distances  in  rough  agreement  with  Hodge  and 
Wallerstein's  value.  Improved  instrumentation  and  more  precise  observations  have  led  to 
the  use  of  many  "secondary"  indicators  in  determining  the  Hyades  distance,  and  some  of 
these  are  explained  briefly  in  the  following  paragraphs. 


Table  2-4.  Hodge  and  Wallerstein  Hyades  distances. 


Indicator 

Distance  (pc) 

Trigonometric  parallaxes 

48.5  ±7.8 

Dynamical  parallaxes  of  visual  binaries 

44.4 

Widths  of  K-emission  lines 

52.1 

Average  of  the  above 

48.3 

Trigonometric  Parallaxes 

The  simplest  of  the  secondary  indicators  is  also  the  most  direct,  as  the  measurement  of 
a star's  trigonometric  parallax  determines  its  distance:  the  inverse  of  that  parallax.  Until 
the  past  two  decades  or  so,  parallaxes  were  of  fairly  low  precision,  with  standard  errors  of 
0".02  or  more.  As  the  distance  to  the  center  of  the  Hyades  lies  somewhere  between  45 
and  50  parsecs  (parallax  between  0".020  and  0".022),  parallaxes  of  Hyades  stars  had  an 
error  roughly  equivalent  to  the  parallax  itself,  sometimes  resulting  in  the  determination  of 
physically  unrealistic  negative  or  very  small  parallaxes.  However,  the  advent  of  automated 
measuring  machines,  the  MAP  (used  at  the  Allegheny  Observatory,  see  Gatewood  et  al., 
1992),  CCD  cameras  and  finally,  the  astrometric  satellite  HIPPARCOS  have  improved  the 
precision  of  modem  parallax  determinations  to  about  0".001,  a gain  of  an  order  of 
magnitude  over  the  past.  This  improved  precision  has  made  possible  fairly  reliable 
estimates  of  the  distance  to  the  cluster  center  from  the  parallaxes  alone,  provided  the 
number  of  stars  considered  in  the  estimate  is  more  than  a handful.  Table  2-5  summarizes 


24 


some  of  the  Hyades  distance  estimates  obtained  by  averaging  the  measured  Hyad 
parallaxes. 


Table  2-5.  Estimates  of  Hyades  distance  incorporating  parallaxes. 


Investigator 

Distance  (pc) 

No.  of  stars 

Upgren  (1974) 

45.7  ±3.8 

11 

Klemola  et  al.  (1975) 

43.5  ± 2.8 

18 

Upgren  et  al.  (1990) 

45.712.1 

23 

Patterson  and  Ianna  (1992) 

45.912.1 

10 

Photometric  Approaches 

Another  class  of  techniques  that  has  enjoyed  a fair  amount  of  popularity  in  recent 
Hyades  distance  determinations  incorporates  photoelectric  photometry  of  Hyades 
members,  usually  intermediate  band  (Mannery  and  Wallerstein,  1971)  or  (R-I)  (Upgren, 
1974;  Eggen,  1982).  Much  has  been  made  of  the  (R-I)  determinations,  as  the  (R-I)  index 
is  not  strongly  dependent  on  metal  abundance.  The  procedure  followed  by  Upgren 
(1974b)  was  to  use  trigonometric  parallaxes  (after  correction  for  known  systematic 
effects)  of  non-member  dwarfs  in  the  Hyades  field  to  establish  a (Mv,  R-I)  main  sequence 


Table  2-6.  Photometric  distance  determinations. 


Investigator 

Colors  used 

Distance  (pc) 

Eggen  (1969) 

(MIf  R-I) 

40.51 1.9 

Mannery  and  Wallerstein  (1971) 

(V,  V-r) 

44.5  1 2.0 

Upgren  (1974b) 

(Mv,  R-I) 

44.0  1 0.8 

Eggen  (1982) 

(Ml  R-I) 

43.6 1 2.0 

Cameron  (1985) 

(V,  B-V) 

44.9  1 2.3 

for  these  background  stars,  which  have  (R-I)  values  between  0.35  and  1.05.  A linear  fit 
was  then  made  to  the  photometry  of  the  Hyades  members  to  establish  a (V,  R-I)  sequence. 
A comparison  of  this  sequence  to  the  (Mv,  R-I)  sequence  for  the  background  stars  yields 
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the  Hyades  distance  modulus,  m-M.5  Eggen  (1969, 1982)  followed  a somewhat  similar 
procedure,  save  that  he  used  (I-Mj)  as  his  distance  modulus,  where  Mj  is  the  absolute  I 

magnitude  of  the  Hyades  stars.  Both  investigators  obtained  a Hyades  distance  in  the 
neighborhood  of  44  parsecs,  although  the  earlier  determination  by  Eggen  gives  a distance 
in  agreement  with  the  ~41  parsec  distance  of  the  classical  convergent  point  solutions. 
There  are  no  obvious  reasons  for  this  discrepancy  between  Eggen's  determinations,  a point 
made  by  Mannery  and  Wallerstein  (1971),  who  obtained  a Hyades  distance  of  44.5  parsecs 
based  on  their  BVr  photometry. 

Cameron  (1985)  came  up  with  a technique  that  used  UBV  photometry  of  Hyades  and 
field  stars  to  obtain  the  cluster  distance  modulus.  After  selecting  a parallax-limited  set  of 
field  stars  redder  than  the  Sun  (0.68  < (B-V)  < 1.00)  from  Gliese's  catalogue,  he  applied 
the  Lutz-Kelker  corrections  and  corrections  for  metallicity,  thereby  obtaining  a zero  age 
main  sequence  (ZAMS).  This  was  compared  to  a Hyades  main  sequence  established  by 
Hyades  dwarfs  with  0.1  < (B-V)  <1.1,  resulting  in  a Hyades  distance  modulus  of  3.26  ± 
0.1 1,  corresponding  to  a distance  of  44.9  ± 2.3  parsecs. 

The  Wilson-Bappu  Effect 

One  of  the  arguments  Hodge  and  Wallerstein  (1966)  presented  against  the  Hyades 
distances  determined  from  the  classical  convergent  point  investigations  involved  the 
Wilson-Bappu  effect.  In  brief,  the  Wilson-Bappu  effect  states  that  there  is  a linear 
relationship  between  the  logarithm  of  the  Call  H and  K emission-line  core  width  and  the 
absolute  visual  magnitude  of  the  star.  Hodge  and  Wallerstein  (1966)  performed  a re- 
determination of  the  parameters  of  this  linear  relationship  using  seven  background  stars 
with  known  parallaxes  and  y Andromedae  A.  Fitting  the  Hyades  stars  to  this  line  resulted 
in  a distance  of  52  parsecs  to  the  cluster  center.  However,  Wilson  (1967)  challenged  this 

5 The  distance  is  obtained  from  the  distance  modulus  by  noting  that  m - M + 5=5  log  d,  where  d is  the 
distance  in  parsecs. 
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value  based  on  his  recalibration  of  the  Wilson-Bappu  effect,  involving  65  stars  with 
"good"  parallaxes,  finding  no  substantial  disagreement  with  a Hyades  distance  of  41 
parsecs.  Lutz  (1970)  settled  the  debate  by  showing  that  Wilson's  linear  regression 
introduced  systematic  errors  into  the  parameters  by  allowing  errors  only  in  the  emission- 
line widths  and  ignoring  those  in  Mv.  The  application  of  an  improved  adjustment 

technique  allowing  errors  in  both  observed  quantities  enabled  Lutz  to  derive  a Hyades 
distance  of  44.3  parsecs,  in  good  agreement  with  other  distance  indicators.  See  table  2-7 
for  a summary  of  these  determinations. 

Table  2-7.  Hyades  distance  estimates  obtained  from  the  Wilson-Bappu  effect. 

Investigator  Distance  (pc) 


Hodge  and  Wallerstein  (1966) 

52 

Wilson  (1967) 

41 

Lutz  (1970) 

44.3 

Techniques  Involving  Binaries 

Using  orbit  solutions  for  several  Hyades  binaries  and  applying  the  local  mass- 
luminosity  relation,  Hodge  and  Wallerstein  (1966)  derived  a Hyades  distance  of  44.4 
parsecs.  McClure  (1982)  followed  a similar  approach,  but  he  was  able  to  directly 
determine  the  mass  of  the  Hyad  vB  22  (HD  27130),  a double-lined  spectroscopic  and 
partially  eclipsing  binary.  After  assuming  a slope  for  the  mass-luminosity  relation,  he  was 
able  to  derive  a cluster  distance  of  49.4  ±1.1  parsecs  from  the  available  data  on  10  visual 
binaries  belonging  to  the  Hyades.  Peterson  and  Solensky  (1987)  revised  this  distance 
down  to  47.0  ±1.2  parsecs,  based  on  a new  estimate  of  the  period  of  one  of  the  binaries, 
Finsen  342,  and  a re-evaluation  of  the  available  data  on  vB  22.  Table  2-8  summarizes 
these  investigations. 
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Table  2-8.  Hyades  distance  determinations  using  binary  star  data. 

Investigator  Distance  (pc) 

Hodge  and  Wallerstein  (1966)  44.4 

McClure  (1982)  49.4  ±1.1 

Peterson  and  Solensky  (1987) 47.0  ±1.2 

Constraints  Imposed  by  Stellar  Models 

With  the  increase  in  computing  power  over  the  past  few  decades,  it  has  become 
possible  to  construct  elaborate  stellar  models  and  produce  theoretical  color-magnitude 
diagrams  and  a mass-luminosity  relation  for  the  Hyades,  which  can  then  be  compared  to 
the  observations.  Problems  became  evident  three  decades  ago  when  Iben  (1963)  and 
Bodenheimer  (1965)  were  unable  to  produce  models  that  matched  the  observed  Hyades 
color-magnitude  diagram  and  the  mass-luminosity  relation  as  derived  from  distances  based 
on  convergent  point  solutions.  Iben  could  achieve  consistency  only  by  assuming  Y=0.26, 
Z=0.0625  in  his  model  calculations,  values  ruled  out  by  spectroscopic  and  photometric 
observations  of  Hyades  members.  The  high  metal  abundance  was  felt  to  be  necessary  to 
explain  the  differences  between  the  observed  Hyades  mass-luminosity  relation  and  that  of 
the  solar  neighborhood.  According  to  Eggen  (1963),  a Hyad  of  0.65  solar  masses  would 
have  a luminosity  equal  to  that  of  the  Sun,  making  the  Hyades  stars  overluminous  for  their 
masses.  Faulkner  (1967)  used  his  quasi- homology  relations  to  point  out  that  an  increase 
in  the  Hyades  distance  from  40.4  parsecs  to  44.4  parsecs  could  resolve  the  differences 
between  the  models  and  observations  without  having  to  postulate  large  metal  abundances. 
The  same  distance  was  obtained  by  Wallerstein  and  Hodge  (1967),  based  on  models 
produced  by  Demarque  (1967).  Later  on,  Hardorp  (1980)  used  his  interior  and 
atmosphere  calculations  to  obtain  a cluster  distance  of  45  parsecs,  and  recently 
VandenBerg  and  Bridges  (1984)  came  up  with  a value  of  49.0  ±1.1  parsecs.  It  should  be 
noted  that  as  the  estimated  Hyades  distance  increased,  the  cluster  mass-luminosity  relation 
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moved  towards  agreement  with  that  of  the  solar  neighborhood.  Table  2-9  is  a short 
summary  of  Hyades  distance  estimates  based  on  stellar  models. 

Table  2-9.  Hyades  distance  constraints  based  on  stellar  models. 


Investigator  Distance  (pc) 

Faulkner  (1967)  > 44. 4 

Wallerstein  and  Hodge  (1967)  > 44.4 

Hardorp  (1980)  45 

VandenBerg  and  Bridges  (1984)  49.0  ±1.1 


Critique  of  Past  Hyades  Distance  Determinations 

As  a glance  at  figure  2-3  will  readily  show,  the  distance  to  the  Hyades  is  still  very 
much  an  open  question,  the  general  trend  being  for  the  cluster  distance  to  increase  with 
time.  Indeed,  if  we  were  to  accept  the  error  bars  in  figure  2-3  at  face  value,  it  would  seem 
that  the  Hyades  cluster  is  receding  from  the  solar  system  at  a speed  of  0.66c,  despite 
having  a measured  radial  velocity  of  ~39  km  s4!  This  is,  of  course,  nonsense,  but  it  serves 
to  illustrate  the  point  that  the  cluster  distance  determinations  are  not  as  accurate  as  the 
given  errors  would  indicate. 

As  the  author's  training  lies  in  the  field  of  astrometry,  it  is  somewhat  difficult  and 
perhaps  inappropriate  to  be  very  critical  of  the  Hyades  distance  estimates  achieved 
through  the  use  of  "secondary"  astrophysical  indicators.  However,  the  fact  that  many  of 
these  indicators  rely  on  trigonometric  parallaxes  in  the  calibrations  is  something  that 
cannot  be  ignored  in  the  evaluation  of  the  precision  or  accuracy  of  their  distance  estimates. 
The  best  currently  available  parallaxes  have  standard  errors  of  0".001,  which  is  more  than 
5%  of  a typical  Hyad  parallax,  while  the  older  parallaxes  used  in  many  calibrations,  such 
as  in  that  for  the  photometric  determinations,  have  quoted  errors  ten  times  this  amount. 
Also,  it  would  seem  that  the  naively  erroneous  application  of  linear  regression  is  not 
restricted  to  just  the  early  convergent  point  investigators.  As  Lutz  (1970)  showed,  Wilson 


29 


(1967)  commited  a similar  mistake  in  his  calibration  of  the  Wilson-Bappu  effect,  thereby 
underestimating  the  distance  to  the  Hyades  by  slightly  over  three  parsecs.  Finally,  it 
should  be  pointed  out  that  some  of  these  indicators,  such  as  early  stellar  models,  use 
observations  of  Hyads  as  a theoretical  check,  making  their  use  as  estimators  of  the  Hyades 
distance  seem  somewhat  like  circular  arguments. 

Given  their  quoted  lesser  precision,  it  is  somewhat  ironic  that  the  distance 
estimates  derived  from  the  secondary  indicators  should  nowadays  be  given  as  much  weight 
as  the  estimates  obtained  from  the  more  precise,  "fundamental"  convergent  point 
solutions.  To  some  extent,  this  has  not  been  due  to  any  deficiencies  in  the  convergent 
point  techniques,  but  is  a consequence  of  the  acknowledgment  that  the  various  sets  of 
proper  motion  and  radial  velocity  data  suffer  from  systematic  errors  that  produce  the  wide 
scatter  observed  in  the  Hyades  distance  estimates.  Nonetheless,  the  convergent  point 
techniques  do  suffer  from  some  serious  deficiencies,  which  are  summarized  in  the 
paragraphs  below. 

All  stars  are  forced  to  have  the  same  velocity.  Even  though  the  velocity  dispersion  of 
the  stars  which  are  members  of  the  Hyades  is  small  (~0.25  km  s1;  see  van  Bueren,  1952, 
and  Gunn  et  al.,  1988),  the  precision  of  modem  radial  velocities  is  approaching  this  value 
(Griffin  et  al.,  1988)  and  shall  soon  reach  a point  where  the  internal  motions  of  the  cluster 
members  will  no  longer  be  "lost  in  the  noise".  As  the  observations  become  precise  enough 
to  discern  the  individual  stellar  motions,  so  too  must  the  data  reduction  techniques  allow 
each  star  to  "go  its  own  way". 

Trigonometric  parallaxes  are  ignored.  Despite  the  gains  in  precision,  trigonometric 
parallaxes  are  left  out  of  convergent  point  solutions  entirely  and  are  relegated  to  the  role 
of  a "secondary"  distance  estimator.  This  should  not  be  the  case,  as  parallaxes  are 
fundamental  observations,  with  just  as  much  importance  as  proper  motions  and  radial 
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velocities.  It  should  be  noted  that  equations  (2-15)  explicitly  contain  the  parallax,6  making 
its  incorporation  into  an  adjustment  procedure  simple  for  someone  familiar  with  the 
appropriate  techniques.  All  that  needs  to  be  done  is  to  use  all  three  of  equations  (2-15), 
rather  than  follow  the  past  procedure  of  ignoring  the  first  and  eliminating  the  parallax  from 
the  second  and  third  equations.  This  way  no  fundamental  information  is  ignored,  and,  as  a 
significant  number  of  Hyads  (>  40)  have  observed  parallaxes,  such  an  adjustment  is  bound 
to  yield  a more  accurate  cluster  distance  estimate. 

Separate  adjustments  are  made  for  each  data  type.  In  the  convergent  point  solutions, 
the  proper  motions  are  used  to  obtain  an  estimate  of  the  convergent  point  location,  and 
then  a separate  adjustment  involving  the  radial  velocities  is  performed  to  obtain  the  cluster 
space  speed.  This  approach  is  outdated;  modem  adjustment  techniques  and  computing 
power  allow  for  an  adjustment  of  all  fundamental  data  simultaneously,  thereby  achieving  a 
more  accurate  estimate.  The  same  criticism  may  be  leveled  at  the  proper  motion  gradient 
methods,  which  incorporate  only  one  proper  motion  component  into  each  adjustment. 

The  data  are  not  used  in  such  a manner  as  to  extract  the  most  information.  One  of  the 
most  serious  criticisms  that  can  be  made  of  the  classical  convergent  point  solutions  is  that 
they  ignore  information  by  considering  only  the  position  angles  of  the  proper  motions  in 
the  estimation  of  the  position  of  the  convergent  point.  The  magnitudes  of  the  proper 
motions  also  carry  information,  and  discarding  them  can  only  result  in  a less  precise 
estimate,  even  though  the  precision  may  seem  better.  This  mistake  is  often  compounded 
by  the  mysterious  and  arcane  ways  in  which  stars  are  selected  for  incorporation  into  the 
adjustment.  For  example,  Schwan  (1991)  presents  145  stars  as  members  of  the  cluster, 
yet  uses  only  62  of  these  to  estimate  the  location  of  the  convergent  point. 

In  light  of  better  observations  and  ever-increasing  computer  power,  it  would  seem 
time  to  retire  the  classical  techniques  in  favor  of  some  other  procedure  that  can 


6 Actually,  the  distance. 
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incorporate  all  fundamental  data  and  relax  the  physically  absurd  requirement  that  all  stars 
must  have  the  same  velocity.  However,  before  these  techniques  can  be  retired,  a better 
procedure  must  be  found  to  take  their  place.  The  purpose  of  this  work  has  been  the 
development  of  just  such  a procedure;  its  outline  and  application  are  presented  in  the 
pages  that  follow. 
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Figure  2-3.  Hyades  distance  determinations 


CHAPTER  3 

LEAST  SQUARES  WITH  PROBABILISTIC  CONSTRAINTS 


Introduction 


The  traditional  approach  in  solving  the  overdetermined1  system  Ma  = x,  where  a is 
the  vector  of  parameters  (unknowns)  and  x is  the  vector  of  observable  quantities 
(observations),  is  through  the  application  of  the  principle  of  least  squares.  In  this 
technique,  we  seek  estimates  of  a such  that 


where  m is  the  number  of  observations  and  n is  the  number  of  parameters.  The  fact  that 
the  magnitude  of  the  innermost  sum  is  squared  gives  the  method  its  sobriquet  of  "least 
squares".  Other  techniques  do  exist,  in  particular  the  method  of  minimum  sum,  which 
obtains  the  parameters  by  setting 


While  the  method  of  minimum  sum  is  more  robust2  than  the  method  of  least  squares, 
making  it  useful  for  identifying  outliers  and  observational  "blunders",  it  is  also  true  that  the 


1 An  overdetermined  system  is  a system  of  independent  equations  in  which  the  number  of  equations 
exceeds  the  number  of  unknowns. 

2 By  "robust",  it  is  meant  that  the  parameter  estimates  obtained  by  the  technique  are  not  strongly 
dependent  on  the  distribution  function  of  the  observations. 


(3-1) 


(3-2) 
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technique  does  not  give  the  most  accurate  estimates  of  a if  x is  free  from  bias3 . 
Furthermore,  the  minimum  sum  does  not  yield  unambiguous  estimates  of  the  parameters 
(Eichhom,1993).  Nonetheless,  Branham  (1990),  who  has  developed  an  algorithm  for  its 
implementation,  suggests  this  technique  because  of  its  usefulness  in  preprocessing. 

In  the  flurry  to  develop  more  robust  adjustment  techniques,  the  method  of  least 
squares  is  occasionally  criticized  as  being  too  specialized,  especially  with  regard  to  the 
distribution  function  of  the  (unknown)  errors.  Many  adjustors  of  data  labor  under  the 
mistaken  impression  that  the  least-squares  technique  requires  the  errors  in  the 
observations  to  obey  the  so-called  normal  law, 


—ft. 


(3-3) 


where  a is  the  standard  deviation  of  the  observations,  and  the  vector  of  (unknowable) 
observation  errors,  is  given  by  ^ = x - x^,  x^  being  the  vector  of  bias-free  observations. 
This  is  not  so,  as  the  principle  of  maximum  likelihood  can  be  invoked  to  show  that  the 
most  accurate  parameter  estimates  are  obtained  by  a least-squares  adjustment  if  the 
distribution  of  errors  is  described  by  any  function  which  shares  its  (unique)  maximum  with 
that  of  the  Gaussian  error  function,  not  just  that  given  by  equation  (3-3).  It  should  also  be 
noted  that  the  Gauss-Markov  theorem,  which  states  that  the  least-squares  estimator  has 
the  minimum  variance  of  the  class  of  linear  unbiased  estimators,  places  no  restriction  on 
the  form  of  the  distribution  function  of  the  errors. 

In  any  adjustment,  whether  by  least  squares  or  some  other  technique,  there  must  exist 
an  adjustment  model,  that  is,  a set  of  equations  relating  the  observations  to  the  unknown 


3 Bias  can  be  defined  mathematically  by  cp(^)  = Ke  20  . where  a is  the  bias.  For  example,  the 

mean,  x , of  a set  of  measurements  of  the  width  of  a box  may  be  thought  of  as  an  unbiased  estimator  of 
this  width,  whereas  x - 5 is  a biased  estimator  of  the  box's  width. 
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parameters.  Often  this  model  is  a simplified  approximation  of  the  actual  relationship, 
which  may  be  unknown,  unknowable  or  deemed  to  be  too  complicated  to  be  considered  in 
the  adjustment.  Whatever  the  reason(s),  any  deficiencies  in  the  model  will  result  in 
"systematic"  errors  in  the  parameter  estimates,  making  it  important  that  the  adjustment 
model  conform  as  closely  as  possible  to  reality.  Many  of  the  "systematic"  differences 
between  observations  obtained  from  different  observatories  or  by  different  instruments  can 
be  explained  by  use  of  inadequate  adjustment  model,  though  it  may  be  impossible  in  many 
instances  to  determine  models  sufficiently  close  to  reality  such  that  the  systematic  effects 
disappear  or  become  unnoticeable. 

General  Least  Squares 


Introduction 

In  general  least  squares,  the  restrictions  imposed  in  the  ordinary  least  squares  (OLS) 
technique  used  by  most  scientists  and  engineers4  are  abandoned.  In  the  OLS  techniques  it 
is  assumed  that 

1)  the  parameters  enter  linearly  into  the  equations  of  condition, 

2)  each  condition  equation  involves  only  one  observation, 

3)  there  exist  no  constraints  among  the  parameters,  and 

4)  the  observations  have  the  same  precision  and  are  uncorrelated. 

Generalized  least-squares  techniques,  such  as  those  outlined  by  Jefferys  (1980)  and 
Eichhom  (1993)  remove  the  above  limitations  by  allowing 

1)  the  equations  of  condition  to  be  nonlinear  with  regard  to  the  parameters. 


4 It  should  be  pointed  out  that  many  scientists  and  engineers  are  unaware  of  the  restricted  circumstances 
required  for  an  OLS  adjustment,  sometimes  resulting  in  the  misapplication  of  the  technique  and 
consequently,  biased  parameter  estimates. 
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2)  each  equation  of  condition  to  contain  more  than  one  observable, 

3)  constraints  (even  probabilistic  ones)  to  exist  among  the  parameters,  and 

4)  the  observables  to  be  correlated  and  of  varying  precision  . 

Note  that  the  relaxation  of  restriction  4)  requires  that  the  covariance  matrix,  o,  of  the 
observations  must  be  known,  at  least  to  within  a constant  factor,  before  the  adjustment  is 
attempted.  This  underscores  the  importance  of  publishing  the  errors  in  the  observations 
along  with  the  observations,  as  the  variances  computed  from  these  errors  comprise  the 
elements  of  a.  Without  these  errors,  it  becomes  practically  impossible  to  achieve  an 
accurate  parameter  estimation  at  first  attempt.  One  can,  of  course,  always  estimate 
variances  and  where  appropriate,  covariances  by  analyzing  the  residuals  from  a preliminary 
adjustment. 

Derivation  of  the  Normal  Equations 


Let  x,)  be  the  vector  of  unbiased  observations,  with  associated  covariance  matrix  a. 
The  observations  are  estimates  of  the  actual,  "true"  values,  x,  and  the  two  are  related 
through  x = x0  + where  £,  is  the  vector  of  corrections  to  the  observations.  We  seek  to 
solve  the  set  of  model  equations,  more  commonly  known  as  the  equations  of  condition, 


F(x0+£,a)  = 0 
G(a)  = 0’ 


(3-4) 


for  the  vector  a of  parameters.  G(a)  is  a vector  of  functions,  which  when  set  to  0, 
constrain  the  adjustment  parameters;  they  are  neither  implicit  nor  explicit  functions  of  any 
observable  quantities  and  establish  the  relationships  (if  any)  among  the  parameters,  a. 
Equations  (3-4)  must  be  solved  for  the  a while  simultaneously  minimizing 


S = 


(3-5) 
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The  solution  is  best  obtained  through  the  application  of  Lagrange  multipliers.  Realizing 
that  the  product  of  F and  G with  any  number  is  0 at  the  solution,  equation  (3-5)  has  as  a 
consequence,  that  also 


S = ^Tcr14-2FrA-2GT5, 


(3-6) 


at  the  solution  (i.e.,  when  F = 0 and  G = 0),  where  A,  S are  yet  undetermined  vectors  of 
Lagrange  multipliers.  The  minimum  of  S is  found  by  setting  to  0 the  partial  derivatives  of 
equation  (3-6)  with  respect  to  the  observations  and  the  parameters.  Thus, 


and 
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or 


and 


2a~%-2XTA  = 0 

-2ArA-2DrS  = 0,  (3-7) 


where  the  Jacobian  matrices  X and  A are  given  by 


X = 


(3-8) 


evaluated  at  the  solutions.  As  o is  nonsingular,  the  first  of  equations  (3-7)  yields  the 
important  result 

5 = oXtA,  (3-9) 

which  is  needed  for  the  computation  of  the  adjustment  corrections  We  use  this 
equation  to  eliminate  £,  from  the  first  of  equations  (3-4),  and  can  now  write  the  set  of 
general  normal  equations  in  the  form 
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ArA  + DrH  = 0 
f(x0  +oXrA,a)  = 0 
G(a)  = 0 


(3-10) 


In  many  cases  of  physical  and  astronomical  interest,  there  exist  no  exact  relations  among 
the  parameters.  In  these  instances,  there  are  no  equations  G(a)  = 0 and  the  normal 
equations  (3-10)  reduce  to 


ArA  = 0 
f(x0  +oXrA,a)  = 0. 


(3-11) 


Equations  (3-10)  and  (3-11)  are  rigorous;  no  approximations  have  been  made.  But,  as 
pointed  out  by  Jefferys  (1980),  they  generally  form  a set  of  unwieldy  nonlinear  equations, 
being  linear  only  when  the  condition  equations  F(x,  a)  = 0 are  linear  with  respect  to  both 
the  observations  and  the  parameters.  Note  that,  as  these  equations  must  be  satisfied  at  the 
solution  (when  a and  ^ have  obtained  their  correct  values),  some  iterative  numerical 
approach  must  be  employed  in  order  to  achieve  a solution. 

Solution  of  the  Normal  Equations 

Both  Jefferys  (1980, 1981)  and  Eichhom  (1993)  present  schemes  for  solving  equations 
(3-11).  These  authors  develop  approaches  based  on  Newton's  method,  though  Jefferys 
(1981)  later  presents  techniques  for  solving  the  normal  equations  via  the  method  of 
steepest  descent,  Marquardf  s algorithm,  or  the  Fletcher-Powell  algorithm.  All  of  these 
methods  require  an  initial  approximation,  a0,  to  the  parameter  vector,  a,  be  available,  with 
Newton's  method  more  often  than  not  requiring  that  the  approximations  be  such  that 


(3-12) 


a 
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in  order  for  convergence  to  be  achieved.  Though  often  readily  available,  there  do  exist 
situations  in  which  good  parameter  approximations  are  nonexistent;  in  these  cases  another 
technique,  e.g.,  one  of  those  outlined  by  Jefferys  (1981)  must  be  used.  Despite  the 
required  closeness  of  the  initial  parameter  approximations,  a procedure  for  solving 
equations  (3-11)  by  Newton’s  method  shall  now  be  presented,  as  it  will  serve  to  illustrate  a 
basic  approach  to  solving  the  nonlinear  normal  equations  and  also  because  it  is  sufficient 
for  our  purposes. 

Let  a0  denote  the  vector  of  approximations  to  the  parameters,  with  a the 
corresponding  vector  of  corrections  to  these  parameters,  i.e.,  a = a0  + a.  In  a similar 
fashion,  Xq  represents  the  vector  of  observations  and  £ is  the  associated  vector  of 
corrections  to  these  observations.  Thus 


Equations  (3-1 1),  developed  as  a Taylor  series  and  broken  off  after  the  first  order  terms, 
then  become 


a = a0  +a 
x = x0+t 


(3-13) 


Ar  A = 0, 


or, 


Fa+X,£  + A,a  = 0 
A0rA  = 0, 


(3-14) 


where  the  subscripts  mean  that  quantities  in  the  condition  equations  and  the  partial 
derivatives  are  evaluated  at  (a0,  Xo).  Note  that  F0  = F(Xq,  a0),  which  should  be  "small"  if 
the  observations  were  precise  and  equations  (3-12)  were  satisfied.  Substitution  of  the 
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second  of  equations  (3-13)  into  equation  (3-9)  yields 


5.  = aXTtA, 


(3-15) 


which  may  then  be  used  in  the  first  of  equations  (3-14),  producing 


F0  + X0oX0  A + A0a  - 0, 


(3-16) 


which  may  now  be  solved  for  A: 


A = - W0  (F0  + A0a), 


(3-17) 


where 


W0=(x0oX0r)_1 


(3-18) 


The  vector  A obtained  from  equation  (3-17)  will  not  be  the  final  A,  as  equation  (3-16)  is 
only  an  approximation.  Call  it  therefore  Ax,  and  find  xv+1  = jq,  + av+1  = av  +av  (v  = 0) 

Now  form  Xv  Aj  and  Fj,  go  through  another  iteration  to  find  A2,  and  0C2,  and  so  forth 

until  further  iterations  no  longer  produce  any  significant  changes.  Of  course, 


where  we  assume  that  the  procedure  has  converged  after  N steps. 

The  elements  of  W,  the  "weight  matrix",  are  well  known  to  anyone  who  frequently 
adjusts  data.  It  is  important  to  note  that  one  is  not  free  to  choose  the  weights  arbitrarily, 
as  has  been  done  in  many  adjustments.  Least  squares  requires  that  the  observational 
weights  be  determined  by  the  Jacobian  matrix  of  the  condition  equations  with  respect  to 
the  observations  and  the  covariance  matrix  of  the  observations;  anything  different  will  lead 
to  biased  parameter  estimates.  In  the  OLS  adjustments  performed  by  most  scientists  and 


N 
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engineers,  X is  an  identity  matrix.  This  results  in  W being  diagonal,  with  the  inverses  of 
the  observational  variances  comprising  the  elements.  Therefore,  the  appropriate  weights 
for  an  OLS  adjustment  are  given  by 


where  o,-  is  the  variance  associated  with  the  i-th  observed  quantity. 

As  with  the  observations,  the  parameter  estimates  are  of  limited  usefulness  unless  they 
are  accompanied  by  their  covariance  matrix.  It  can  be  shown  that  this  matrix  is  given  by 


The  reader  can  find  a derivation  of  equation  (3-20)  in  most  advanced  statistics  texts. 
Branham  (1990)  gives  a simple  development  for  linear  systems,  whereas  Duane  Brown 
(1955)  and  Jefferys  (1980,  1981)  present  a more  general  treatment.  The  standard  errors 
of  the  parameters  are  equal  to  the  square  roots  of  the  diagonal  elements  of  this  matrix,  or 


A Simple  Example 

The  matrix  notation  used  in  presenting  data  adjustment  schemes  is  similar  to  the 
terminology  employed  by  astrometrists;  both  seem  to  contribute  to  the  confusion  of  the 
novice,  even  though  the  concepts  are  actually  quite  simple.  It  would  therefore  seem 
appropriate  to  illustrate  the  general  least-squares  procedure  outlined  in  the  previous 
section  by  use  of  a physical  example,  one  that  has  somewhat  of  a bearing  on  the  Hyades' 


(3-19) 


(3-20) 


(3-21) 


distance. 
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As  mentioned  in  chapter  2,  the  Wilson-Bappu  effect  states  that  there  is  a linear 
relationship  between  the  logarithm  of  the  Call  emission  line  width  and  the  absolute  visual 
magnitude  of  the  star.  We  may  therefore  write 

My  -q\ogW  + b = 0,  (3-22) 


where  Mv  is  the  star's  indirectly  observed  absolute  visual  magnitude,  W is  the  observed 
width  of  the  emission  line,  and  q and  b are  the  unknown  parameters.  The  equation  above, 
when  set  up  for  each  of  the  stars  involved  in  the  adjustment,  form  the  set  of  condition 
equations 


^Mv  -q\ogWx  +b  ' 


MVj  -q\ogW2  +b 


= 0, 


(3-23) 


MVm-q\ogWm+b^ 


where  the  subscripts  indicate  the  particular  star  considered,  m being  the  total  number  of 
stars.  These  quantities  are  considered  to  be  the  observations,  resulting  in 

*L=(m„  Mv,  W2  - M,_  W.)  (3-24) 


Note  that  each  star  generates  two  "observations"  and  that  all  these  are  incorporated  into 
one  vector.  Accompanying  each  Mv  and  W should  be  the  associated  standard  error, 
denoted  by  eM  and  ew,  respectively.  Assuming  the  observations  to  be  uncorrelated,  let  the 
covariance  matrix  of  these  observations  be  given  by 


e2Mi  0 0 0 " 

0 4,  0 0 

0 0 £2Mb  0 

0 0 • • • 0 EWm  ^ 


(3-25) 
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We  now  set  up  the  Jacobian  matrices  A and  X by  taking  the  partial  derivatives  of 
equations  (3-28)  with  respect  to  the  parameters  and  the  observations.  This  results  in 


'-logw; 

-logW2 


f 

1 


and 


l-logWm  1) 


gln(10) 


0 

<y  ln(10) 


0 0 
0 0 


0 0 0 0 

V 


gln(10) 


as  the  parameter  vector,  a0,  is  given  by 

< ={q  b) 


(3-26) 


(3-27) 


(3-28) 


Having  constructed  the  matrices  needed  to  perform  the  adjustment,  the  parameter 
estimates  can  be  obtained  by  proceeding  as  follows: 

1)  Start  with  v = 0,  i.e.,  and  assume  a plausible  a£  = (q0  b0),  where  q0  and  b0  are 
initial  guesses  for  the  parameters,  obtained  from  Eggen  (1967).  Set  up  the  diagonal 
covariance  matrix  a from  the  observational  errors. 

2)  Compute  A^,  Xv,  Wv,  and  Fv  from  the  current  values  of  jq,  and  av. 

3)  Calculate  the  vector  of  the  parameter  corrections,  cq,  from 

a = -(A^WvAv)_1A^WvFv 
Av=-Wv(Fv  + Avav) 

= <AV 
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4)  Compute  the  new  quantities  x^,  av+1  and  Fy+1 

5)  If  the  components  of  Fv+1  are  sufficiently  small,  go  to  step  6)  below.  Otherwise, 
go  to  step  2). 

6)  Convergence  has  been  achieved;  compute  the  errors  in  the  parameters  by  using 


to  implement,  requiring  slightly  more  work  than  that  involved  in  an  ordinary  least-squares 
approach.  Conceptually,  this  is  indeed  true;  however,  the  implementation  of  a general 
least-squares  adjustment  requires  considerably  more  planning  and  far  more  computer 
power  than  the  corresponding  OLS  technique.  The  reason  for  this  lies  in  the  large  size  of 
the  matrices  associated  with  the  observations.  In  the  simple  example  presented,  X is  of 
the  order  of  m x 2m,  whereas  A is  only  m x 2 in  size.  An  adjustment  involving  a more 
elaborate  set  of  condition  equations  and  many  observations  can  easily  produce  a matrix  X 
with  over  a million  elements,  making  XgXt  of  similar  size.  This  product  must  be  inverted 
in  order  to  obtain  the  "weight"  matrix,  W,  a herculean  task  for  even  modem  computers. 
Fortunately,  X is  "sparse",  that  is,  most  of  its  elements  are  0.  This  results  in  XoXT  being 
block  diagonal  and  symmetrical,  enabling  the  computation  of  W by  the  inversion  of  the 
much  smaller  individual  blocks  via  efficient  inversion  routines.  There  also  now  exist 
matrix  routines  (e.g.,  Branham,  1990)  that  can  greatly  facilitate  multiplications  and  other 
operations  involving  sparse  matrices,  and  these  should  be  seriously  considered  for 
inclusion  in  the  adjustment  program. 


<*„  =(a'wa)'' 


It  would  seem  from  the  above  that  a general  least-squares  solution  is  not  at  all  difficult 
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Probabilistic  Constraints 

Eichhom  (1978)  has  further  extended  the  theory  of  least  squares  to  enable  the 
estimation  of  parameters  which  are  samples  of  a Gaussian  distribution  with  known 
covariance  but  which  whose  modes  (i.e.,  average  values)  are  unknowns.  There  are 
numerous  astrometric  instances  where  such  a technique  would  be  valuable,  one  example 
being  the  restriction  of  certain  parameters  in  an  overlap  solution  incorporating 
photographic  plates  taken  with  the  same  instrument.  In  this  instance,  the  focal  length  of 
the  optical  system  can  be  allowed  to  vary  about  some  initially  undetermined  mean  value, 
resulting  in  more  accurate  star  positions  than  would  be  the  case  if  the  focal  length  were 
unrestricted  (a  "free"  parameter)  or  forced  to  have  the  same  value  for  all  plates.  Another 
example  involves  restricting  the  velocities  and  parallaxes  (distances)  of  the  member  stars 
of  an  open  cluster  to  obey  Gaussian  distributions  with  known  variances,  and  is  the  focus 
of  this  work. 

Least  squares  with  probabilistic  constraints  assumes  the  presence  of  two  parameter 
types,  the  first  of  which  is  that  vector  of  unconstrained,  "conventional"  parameters  familiar 
to  least-squares  practitioners.  The  second  type  of  parameter  is  the  vector  of  parameters 
subject  to  probabilistic  constraints;  hence  the  technique's  name,  "least-squares  with 
probabilistic  constraints".  This  vector  of  constrained  parameters  really  consists  of  two 
parts;  the  vector  g of  the  mode  and  the  random  vector  P of  the  individual  variations.  We 
assume  that  we  have  an  approximation  gj  to  gT  =(gl  g2  • ■ • gk)  and  thus  need  to 

find  the  k- vector  y of  corrections,  implying  that  we  have  put 

g = g0+Y-  (3‘29) 

In  addition  to  each  component  gK  of  g belongs  a |iK- vector  P*  = (PKl  P^  • • • Pk^  ) > 
representing  the  |iK  randomly  distributed  samples  PK  of  the  component  gK  of  the  vector  g. 

k 

We  combine  all  the  (iK- vectors  PK  into  the  Z|1K -vector  P,  i.e., 

K=1 
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pt = (pr  pi  - pi) 


and  introduce  the  £|lK -vector  b by 


K=1 


br=(ei£i+Pi  •••  ekg*+pk), 


(3-30) 


where  the  pK- vector  eK  is  given  by 


=(i  - i). 


We  assume  that  the  £|XK  components  of  P obey  the  distribution  function  (Eichhom, 

K=1 


1978), 


_JLprp->p 

¥(p)  = & 2 


(3-31) 


or  some  other  exponential  function  that  shares  the  same  maximum. 
Let  the  equations  of  condition  be  given  by 

F[x,  b(g,  P),  a]  = 0 


(3-32) 


where,  as  before,  a = a0  + a.  The  joint  distribution  function  of  the  observations  and  the 
restricted  parameters  is  then 


Y(P)<P(^)  = ^ 2 e2  -Ke 


(3-33) 


We  wish  to  solve  equations  (3-32)  for  a and  g while  simultaneously  minimizing 


5 = ^a-^  + prp-'p 


(3-34) 


5 Note  that  here  we  have  "sneaked"  in  the  assumption  that  both  distribution  functions  are  Gaussian. 
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As  before,  the  technique  of  Lagrange  multipliers  may  be  used,  resulting  in 

S'  = + PVP  - 2Ar(X^  + Aa  + BP  + Cy  + F0 ),  (3-35) 


where  B is  the  Jacobian  matrix  of  partial  derivatives  of  F with  respect  to  the 
probabilistically  restricted  parameter  corrections,  P,  and  C is  the  matrix  of  partial 
derivatives  of  F with  respect  to  the  means  of  the  restricted  parameters,  g,  i.e.. 


B = 


'dF\ 

UgJ 


(3-36) 


Note  that  the  matrices  B and  C are  built  up  of  the  same  components  but  that  the 

k 

dimension  of  C is  m x k,  while  that  of  B is  m x X|lK.  Taking  the  partial  derivatives  of  the 

K=1 

above  and  setting  them  equal  to  0 gives  the  normal  equations 


ArA  = 0 

ArC  = 0 (3-37) 

F(x,a,g)  = 0, 


in  which  and  P have  been  expressed  in  terms  of  A by  the 

£ = oXr  A 
P = pBrA 


(3-38) 


Solving  the  system  for  A yields 

A = -W(A0a  + C0Y  + F0),  (3-39) 

where  the  weight  matrix  W is  now  given  by 


W = (xoaXj  + BoPB  J )”' 


(3-40) 
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Equation  (3-40)  is  a beautiful  illustration  of  the  similarity  between  the  observations  and 
the  restricted  parameters.  Both  follow  Gaussian  distributions,  both  their  covariance 
matrices  must  be  known,  and  both  enter  into  the  mathematics  of  the  adjustment  in  like 
fashion. 

The  parameter  corrections,  a and  y,  may  be  obtained  by  solving  the  system 


'a;wa, 

at0wcA 

f°V 

( S.T\ 

WF0  =0 

(3-41) 

,c;wa, 

cjwc0J 

IyJ 

CT 

As  the  normal  equations  will  be  generally  nonlinear,  a numerical  approach  similar  to  that 
outlined  in  the  previous  section  must  be  used  to  obtain  a solution,  a,  g and  x are  updated 
according  to 

a = a0  +a 

g = g„+Y  (3‘42) 

x = x0  + £, 


where  E,  is  given  by 

^ = -oXr  W(A0a  + C0y  + F0)  (3-43) 

Once  convergence  has  been  achieved,  the  vector  |J  is  calculated  from 

P = ~pBT W(  A0a  + C0y  + F0 ) , (3-44) 

and  then  the  restricted  parameters  from  equation  (3-30).  Finally,  the  covariance  matrix  of 
the  parameters  is  obtained  from  the  analogy  to  equation  (3-20). 
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The  Condition  Number  and  the  Efficiency 

Occasionally  a practitioner  of  least  squares  will  hear  some  matrix  described  as  "ill- 
conditioned".  This  has  nothing  to  do  with  the  everyday  conventional  definition  of 
"condition",  as  a matrix  has  no  health;  rather,  it  refers  to  the  suitability  of  the  matrix  for 
inversion.  This  suitability  is  measured  by  the  "condition  number,"  k,  of  the  matrix,  defined 
for  symmetric  positive  definite  matrices  by 


where  ^min  and  X,max  are  the  minimum  and  maximum  eigenvalues  of  the  matrix.  The  larger 
the  condition  number,  the  less  accurate  the  inverse.  This  is  because  k is  essentially  the 
decadic  logarithm  of  the  factor  by  which  the  inaccuracies  in  the  elements  of  the  matrix  are 
multiplied  when  calculating  the  inverse  (see  Branham,  1990  for  more  details).  For 
example,  a k of  the  order  of  106  would  indicate  that  a matrix  could  not  be  inverted  using 
single  precision  arithmetic,  which  has  a precision  of  lO6. 

The  efficiency,  e,  of  a symmetric,  positive  definite  matrix,  Q,  defined  by  Eichhom 


where  n is  the  rank  of  Q,  measures  still  another  property  of  the  matrix.  If  Q is  the 
covariance  matrix  of  a vector  a,  the  efficiency  is  a measure  of  the  degree  of  correlation 
between  the  elements  of  a.  A low  efficiency  indicates  a high  degree  of  correlation, 
originating  from  the  off-diagonal  elements  of  Q,  whereas  a high  efficiency  (e  ~ 1)  indicates 
a low  degree  of  correlation  between  the  elements  of  a.  The  application  of  this  to  least- 
squares  adjustments  should  be  obvious,  as  the  efficiency  of  the  covariance  matrix  of  the 


(3-45) 


(1993)  as 


(3-50) 
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parameters  will  indicate  the  degree  of  correlation  among  the  parameters.  It  is  important  to 
note  that  the  efficiency  and  the  condition  number  are  independent  of  each  other. 


CHAPTER  4 

THE  HYADES  DISTANCE 


The  Normal  Equations  for  the  Hvades 


Having  outlined  the  theory  behind  least  squares  adjustments  with  probabilistic 
constraints  in  chapter  3,  it  is  now  time  to  focus  on  the  application  of  this  theory  to  the 
problem  of  the  determination  of  the  distances  and  velocities  of  the  stars  belonging  to  the 
Hyades  cluster1 . The  first  step  towards  this  problem's  solution  is  the  same  as  the  initial 
step  in  any  adjustment  procedure,  namely,  to  obtain  the  observations  and  set  up  the 
equations  of  condition.  The  latter  has  already  been  done,  for  in  chapter  2 we  saw  that  for 
any  star, 


f . \ 

nr 

acos8 

8 


= 7tR,(-5)R2(a)x,  or 


f r-xcosacos8-ysinacos8-zsin8 
iqia  + xjtsina-  yttcosa 
jcp5  + xTt  cosa  sin  5 + yit  sin  a sin  8 - zjt  cos8 


= 0, 


(4-1) 


where  we  have  replaced  the  proper  motion  components  a cos  8 and  8 by  |ia  and  }i5,  the 
units  being  the  only  difference.  As  before,  r is  the  star's  radial  velocity,  a,  8 are  the 
components  of  its  position,  7t  is  its  parallax  (in  " A.U.1),  and  x,  y,  and  z are  the  stellar 
equatorial  velocity  components  in  the  same  units  as  r . The  set  of  equations  (4-1) 


1 Indeed,  in  the  first  paper  published  on  least  squares  with  probabilistic  constraints,  Eichhom  (1978) 
suggests  the  application  of  this  technique  to  open  clusters. 
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involving  all  stars  used  in  the  adjustment  shall  form  the  condition  equations,  F.  Therefore, 


r,  -xcosa,  cos8j  -ysmoq  cos8!  - zsm8, 

K|ial  +X7tj  sina,  — yTCj  cosoq 
Kja.si  H-xTij  cosa,  sinS!  + ynl  sinoq  sinSj  - zjtj  cos5j 

rm  - xcosam  cos5m  - y sinam  cos8m  - zsin8m 
W™  + xttm  sinam  - yttm  cosam 
JVs*.  + cosam  sin8„  + yttm  sinam  sin5m  - zttm  cos5m , 


where  m is  the  number  of  stars. 

Obviously,  the  observations  for  each  star  are  its  proper  motion  components,  radial 
velocity  and  parallax.  In  most  other  astrometric  circumstances,  the  position  (a,  5)  is  also 
regarded  as  a pair  of  observations;  however,  in  this  case  the  positions  are  known  with  such 
great  precision  (compared  to  the  errors  in  the  other  observables)  that  they  may  treated  as 
known  parameters  and  not  subjected  to  adjustment.  It  should  be  pointed  out  that,  though 

the  Hyades  have  been  an  object  of  study  for  many  decades,  a full  set  of  all  four 
observables  xj  = (q  |i(I.  p.g.  7t),  is  available  for  only  a small  number  of  the  stars  in  the 

cluster.  Fortunately,  the  adjustment  procedure  can  accommodate  even  the  minimum  set 
xo  = (M-o^  (%)>  by  simply  leaving  out  those  components  of  F which  involve  the  star's  radial 

velocity,  i.e.,  the  first  of  equations  (4-2).  In  the  discussion  to  follow,  it  will  simplify 
matters  to  assume  that  all  four  observables  are  available  for  each  star,  so  that  the  vector  of 
observations  is 


(4-3) 


with  associated  covariance  matrix,  a: 


CJ= 


(4-4) 
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0 0 


V 


0 0 
0 0 


where  the  e denote  standard  errors  in  the  appropriate  quantities.  The  above  assumes  that 
the  observations  are  uncorrelated;  correlated  observations  will  have  non-zero  off-diagonal 
elements.  Note  that  the  observation  vector  is  of  length  4m  and  that  the  sparse  observation 
covariance  matrix  has  dimensions  4m  x 4m. 

The  unknown  parameters  in  equations  (4-2)  are  the  three  cartesian  velocity 
components  and  the  parallaxes  of  the  stars,  the  parallaxes  entering  into  the  adjustment  as 
both  observations  and  parameters.  All  parameters  shall  be  stochastically  (probabilistically) 
constrained,  the  reasons  for  their  restriction  being  set  forth  by  Eichhom: 

It  is  inaccurate  to  assume  that  all  stars  have  the  same  velocity.  Their  velocities  will 
rather  obey  a certain  distribution  function,  for  which  an  (if  necessary,  trivariate) 
Maxwellian  distribution  function  would  be  a good  approximation. . . . The  estimation 
of  the  parallaxes  of  the  member  stars  of  a moving  cluster  poses  a similar  problem.  If 
these  are  estimated  individually  without  any  restrictions,  they  will  have  more  or  less 
large  errors  and  thus  the  cluster  will  seem  more  elongated  in  the  line  of  sight  than  it 
actually  is.  (Eichhom  , 1978,  p.  355-356)2 

Therefore  the  vector  of  conventional  unrestricted  parameters,  a,  does  not  exist  in  our 
problem,  and  we  may  write 


2 Note  that  in  fact,  the  velocity  distribution  function  cannot  be  Maxwellian,  because  the  velocity  of  a star 
at  the  cluster's  boundary  cannot  have  a radial  component  (with  respect  to  the  cluster  center). 

Nevertheless,  it  is  better  to  assume  a Maxwellian  distribution  than  to  constrain  the  velocity  to  be  the  same 
for  each  star. 


gr  = ((x)  (y)  (z)  <tc)  •••  (x)  (y)  (z)  (jc)) 

+ (Pl  ^2  P3  P4  P4m-3  p4m-2  ^4^-1  P4m)» 


(4-5) 


54 


where  the  P's  are  the  stochastic  parameters,  and  the  expectation  values  of  the  parameters 
are  given  by 

(x)  = x0  + yp  (y)=y0+Y2,  <z}  = z0 +y3,  (rc>  = jc0 +y4> 

the  Ys  being  the  parameter  corrections.  The  expectation  values  yield  the  average  cluster 
velocity  and  parallax;  these  are  added  to  the  stochastic  P's  (which  have  expectation  value 
zero)  to  yield  the  parameters  for  each  star,  y is  a 4- vector,  whereas  P is  a 4m  vector; 
hence  the  form  of  equation  (4-5).  Before  attempting  the  adjustment,  reasonable  values  for 
the  elements  of  the  covariance  matrix,  p,  of  the  constrained  parameters  must  be 
determined  by  some  means.  As  the  parameters  are  assumed  to  be  uncorrelated,  this  matrix 
is  of  the  form 

ef*  0 •••  0 0^ 

0 efj,  0 0 

p=  : , (4-6) 

0 0 eL  0 

,0  0 •••  0 

where  the  e's  refer  to  the  dispersions  in  the  indicated  parameters.  It  should  be  pointed  out 
that  the  dispersions  in  the  cartesian  velocity  components  can  be  assumed  equal,  and  that 
the  value  of  this  dispersion  can,  at  least  in  principle,  be  obtained  from  a study  of  the  radial 
velocities  of  the  stars  which  are  Hyades  members.  The  dispersion  in  the  parallaxes  can  be 
obtained  from  the  angular  size  of  the  cluster  or  from  dynamical  theory. 

We  now  set  up  the  various  Jacobian  matrices  required  for  the  adjustment.  An 
application  of  equations  (3-8)  and  (3-40)  to  the  equations  of  condition  results  in  the 
matrices 


^-cosotj  cos  8, 

- sin  cos  8, 
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- sin  8, 

0 

Ttsina, 

-jt  cos  a, 

0 

x sin  a, -y  cos  a, 

c = 

xccosa,  sin  8, 

7t  sin  a,  sin  8! 

-TtCOsSj 

xcosa,  sin8j  + y sinctj  sin  8;  -ZCOS8, 

jtcosam  sin  8„ 

7tsina„  sin  8,, 

-7CCOsSm 

x cos  a„  sin  8„  + y sin  am  sin  8m  — z cosS, 

T 0 0 0 

0 10  xsino^-ycosa, 

0 0 1 xcosa^inSj  + ysina^inS^zcosS, 


0 0 0 0 
0 0 0 0 
^0  0 0 0 


0 

0 

0 

xsina„-ycosa„ 

x cosa„  sinS,,  + y sina„  sin8m  - z cos8m 


Note  that  X is  sparse  and  block  diagonal,  each  block  being  3 x 4 in  the  case  of  a full 
complement  of  observations  for  the  star,  giving  X a maximum  size  of  3m  x 4m.  If  the 
parallax  is  not  available,  the  block  reduces  to  a 3x3  identity  matrix,  and,  if  only  proper 
motions  have  been  observed  for  the  star,  the  block  is  smaller  still  (a  2 x 2 identity  matrix). 
The  lack  of  an  observed  parallax  has  no  influence  on  C (which  has  maximum  dimensions 
of  3m  x 4),  though  a star  with  only  proper  motions  available  will  have  just  two  of  the  three 
rows,  as  the  first  of  equations  (4-1)  must  be  dropped.  Equations  (3-40)  also  establish  the 
matrix  B as 


C.  = 


-cosa,  cos8;  -sina^cosS, 
7tsma,  -7t  cos  (X; 
k cosa,  sin  8,  n sin  a,  sin  8. 


-sin  8,  0 

0 xsina, -ycosot; 

— 7t  cos8,  x cosa,  sin  8,  + y sin  a,  sin  8;  - z cos8, 


Cj  0 0 > 

0 C2  0 


{0  0 ...  CJ 


(4-8) 


There  being  no  conventional,  unconstrained  parameters,  the  solved  normal  equations 


(3-41)  reduce  to 


(4-9) 
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Y = -(crWc)"1CrWF0) 

where  W and  F0  are  defined  by  equations  (3-18)  and  (3-10).  An  iterative  procedure 
similar  to  that  outlined  in  chapter  3 is  then  used  to  obtain  the  vector  y.  Once  convergence 
has  been  reached,  the  parameter  covariance  matrix  is  found  from  equation  (3-25)  and  the 
vector  p calculated  from  equation  (3-44).  Finally,  the  stochastically  restricted  parameters 
(the  velocity  components  and  parallaxes  of  each  star)  are  computed  from  equation  (4-5). 

Application  of  the  Technique  to  Model  Clusters 
Description  of  the  Models 

One  of  the  truly  wonderful  aspects  of  modem  electronic  computers  is  their  capability 
to  realistically  simulate  many  physical  situations,  including  the  motions  of  stars  which  are 
members  of  open  clusters.  This  makes  it  possible  to  construct  a hypothetical  cluster  and 
generate  "observations"  of  its  member  stars,  which  may  then  be  used  in  various  adjustment 
procedures  to  test  the  accuracies  of  their  parameter  estimates.  Nature  does  not  afford  us 
this  luxury,  as  astrometrists  have  no  way  of  knowing  the  true  velocities  and  locations  of 
stars  in  an  open  cluster.  The  introduction  of  a new  technique  naturally  raises  questions  as 
to  its  ability  to  determine  more  accurate  parameter  estimates  than  those  obtained  by 
previous  methods,  questions  that  cannot  be  easily  resolved  by  application  to  the  real 
world.  Therefore,  it  would  seem  advisable  to  utilize  some  of  the  available  computing 
power  to  generate  model  clusters  to  serve  as  tests  for  the  method  of  least  squares  with 
probabilistic  constraints. 

There  is  no  need  to  perform  lengthy  n-body  calculations  to  generate  the  models,  as  the 
important  consideration  is  that  the  velocity  dispersion  of  the  "stars"  should  be  small 
compared  to  the  average  cluster  space  velocity.  The  adopted  equations  of  condition  do 
not  require  that  a star  be  gravitationally  bound  to  the  other  stars  in  the  cluster.  This 
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simplification  makes  for  rapid  generation  of  models  with  a large  number  of  stars,  even  on 
"slow"  computers. 

Twenty-four  model  "associations"  of  100  stars  each  were  constructed  and  placed  into 
three  distinct  sets.  The  first  set  consisted  of  ten  models  having  centers  located  at  a 
distance  of  50  parsecs  (parallax  of  0".02)  at  equatorial  position  (60°,  15°).  The  distances 
of  the  individual  stars  from  the  cluster  centers  were  allowed  to  vary  according  to 

rc  = 5 md(0,  1), 

where  rc  is  the  distance  in  parsecs  of  the  star  from  the  cluster  center,  and  md  (0,  1) 
represents  a Gaussian  random  number  between  0 and  1.  A distribution  of  distances  from 
the  cluster  center  belonging  to  one  of  the  models  is  displayed  in  figure  4-1,  and  figure  4-2 
is  a plot  of  the  stars'  equatorial  positions. 

The  average  space  velocity  for  these  ten  association  models  has  the  components 
vT  = (-1.4,  44.5,  6.5)  km  s4  and  thus  a speed  of  45  km  s4.  This  velocity  vector  yields  a 
convergent  "point"  position  at  (91°. 8,  8°.3).  As  with  the  distances,  the  velocity 
components  of  the  individual  stars  were  allowed  to  deviate  from  the  average  values 
according  to  a Gaussian  distribution,  the  differences  being  given  by 

v =0.5  md  (-1,1) 

Vyc  = 0. 5 rnd  (-1,1), 
v2e  =0.5  rnd  (-1,1) 

where  md  (-1,1)  represents  a Gaussian  random  number  between  -1  and  1. 

Next,  the  "observations"  were  calculated  by  means  of  equations  (2-14).  These 
observations  were  then  "dirtied  up"  by  the  addition  of  synthetic  errors  produced  by 
Gaussian  random  number  generators.  The  width  of  the  Gaussian  used  to  generate  the 
proper  motion  errors  was  set  at  0.005  " yr4,  with  corresponding  values  of  0.5  km  s4  and 
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Distance  from  cluster  center  (pc) 

Figure  4-1.  Distribution  of  stellar  distances  from  cluster  center. 
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0".01  for  the  radial  velocity  and  parallax  errors,  respectively.  The  effects  of  these  errors 
upon  the  proper  motions  and  radial  velocities  is  minor;  this  is  not  the  case  with  the 
parallax  "observations",  as  the  standard  deviation  of  0".01  is  one-half  the  value  of  the 
mean  parallax.  Figure  4-3  depicts  the  number  distribution  of  the  parallax  "observations" 
after  the  errors  have  been  applied;  note  the  huge  scatter  in  the  parallaxes,  some  of  which 
have  even  become  negative. 


Parallax  (arcsec.) 

Figure  4-3.  Frequency  distribution  of  synthetic  parallax  "observations". 

In  order  to  test  the  sensitivity  of  the  adjustment  procedure  to  the  distance  and  velocity 
distributions  of  the  stars  within  a cluster,  the  second  group  of  ten  model  associations  was 
generated  with  the  same  center  location  and  average  velocity  as  those  of  the  first  set. 
However,  the  distances  of  the  stars  from  the  cluster  center  and  the  velocities  were  chosen 
at  random  from  a uniform  distribution,  rather  than  from  a Gaussian.  In  order  to  keep 
consistency  with  the  previous  group,  the  maximum  distance  from  the  cluster  center  was 
set  at  5 parsecs  and  the  velocity  components  were  allowed  to  deviate  from  the  average 
value  by  no  more  than  0.5  km  s'1.  Synthetic  "observations"  were  then  produced  and 
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Gaussian  errors  applied  in  the  fashion  described  for  the  first  set  of  associations.  Figure  4- 
4 is  a depiction  of  the  equatorial  positions  of  the  stars  belonging  to  one  model  in  this 
second  group  of  associations;  note  that  the  cluster  is  more  compact  than  its  counterparts 
from  the  Gaussian  set. 


The  third  set  of  models  consisted  of  four  associations  constructed  in  a fashion  identical 
to  those  in  the  first  group,  save  that  the  position  of  the  cluster  centers  was  moved  from 
(60°,  15°)  to  (110°,  45°).  These  associations  also  had  their  average  velocity  altered  to 
vT  = (-38.4,  6.8,  22.5)  km  s4,  yielding  a convergent  "point"  at  equatorial  position  (170°, 
30°)  and  a \ of  489.7  for  the  cluster  center.  The  primary  purpose  for  the  creation  of  this 
smaller  third  group  was  to  make  sure  that  the  adjustment  techniques  were  not  sensitive  to 
the  direction  of  the  cluster's  velocity  (and,  consequently,  to  the  position  of  the  convergent 
point)  and  to  its  position  in  the  sky. 
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Evaluation  of  the  Adjustment  Procedures 

Using  the  equations  and  procedures  outlined  in  chapter  3 and  in  the  first  part  of  this 
chapter,  least-squares  adjustments  involving  probabilistic  constraints  were  performed  on 
all  model  associations,  with  the  parallax  dispersion  set  at  0".006  and  the  dispersions  in  the 
velocity  components  set  equal  to  0.5  kms1.  The  global  parameter  estimates  for  the 
Gaussian  association  models  are  summarized  in  table  4-1.  It  is  readily  apparent  that  there 
is  some  scatter  in  the  average  association  distance  and  velocity;  however,  the  averages  of 
these  estimates  are  comfortably  close  to  their  true  values.  Also  note  that  no  model  has 
global  parameters  that  deviate  from  the  true  values  by  more  than  the  quoted  standard 
error. 

Table  4-1  is  reassuring  in  that  it  presents  evidence  that  an  adjustment  involving 
stochastic  constraints  can  reliably  characterize  the  gross  characteristics  of  an  association 
or  cluster.  However,  the  thrust  of  the  stochastic  constraint  technique  is  that  the  velocities 
and  parallaxes  of  the  individual  stars  can  be  directly  determined  without  the  errors  in  the 
observed  quantities  exercising  their  full  distorting  influence,  and  it  is  essential  to  see  how 
closely  these  estimates  conform  to  reality.  Some  indication  can  be  gleaned  from  table  4-2, 
which  gives  the  average  of  the  absolute  deviations  of  the  stochastic  velocities  and 
parallaxes  from  the  stars'  true  values.  It  would  seem  that  we  can  match  a star's  velocity  to 
within  1 km  s1,  and  the  parallax  to  about  0".001,  which  is  the  quoted  precision  of  the 
MAP  parallaxes  (Gatewood  et  al.,  1992).  For  a mean  distance  of  50  parsecs,  an  error  of 
0".001  in  the  parallax  would  give  a corresponding  distance  error  of  about  2.5  parsecs. 
This  does  not  seem  very  impressive;  however,  as  stated  previously,  these  adjustments 
were  performed  to  compare  to  the  results  of  the  adjustments  involving  the  uniformly- 
distributed  model  associations.  Therefore,  the  constraint  in  the  parallax  was  deliberately 
set  to  a very  loose  value  of  0".006  (15  parsecs).  It  shall  be  seen  later  that  an  appropriate 
"tightening"  of  the  parallax  restriction  can  result  in  an  average  absolute 
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Table  4-1.  Global  parameters  for  Gaussian  association  models 


Model  No. 

vT„  (km  s*9 

v„„  (km  S' 9 

v-  (km  S' 9 

Kr  (arcsec.) 

1 

-1.29  ± 1.17 

44.58  ± 0.76 

6.47  ± 0.25 

0.0199  ± 0.0017 

2 

-0.41  ± 1.18 

43.96  ± 0.77 

6.48  ± 0.27 

0.0209  ± 0.0018 

3 

0.36  ± 1.09 

43.49  ± 0.69 

6.7  ± 0.22 

0.0217  ± 0.0017 

4 

-1.85  ± 1.2 

44.94  ± 0.78 

6.39  ± 0.24 

0.0195  ±0.0016 

5 

-2.71  ± 1.35 

45.42  ± 0.85 

6.34  ± 0.28 

0.0186  ±0.0017 

6 

-0.93  ± 1.13 

44.29  ± 0.74 

6.56  ± 0.23 

0.0203  ± 0.0017 

7 

-2.96  ± 1.26 

45.54  ±0.81 

6.31  ±0.25 

0.0183  ±0.0016 

8 

-0.9  ± 1.13 

44.3  ± 0.73 

6.43  ± 0.25 

0.0204  ± 0.0017 

9 

-2.38  ± 1.16 

45.25  ± 0.74 

6.49  ± 0.24 

0.0188  ±0.0016 

10 

-0.74  ± 1.18 

44.16  ±0.74 

6.5  ± 0.22 

0.0207  ± 0.0017 

Average 

-1.38  ±1.07 

44.59  ±0.68 

6.47  ±0.11 

0.0199  ±0.0011 

Table  4-2.  Deviations  of  stochastic  parameters  from  true  values  for  Gaussian  models. 
Model  No.  1 Av,  1 (km  s1)  1 Avv  1 (km  s1)  1 Av,  1 (km  s1)  1 Ak  1 (mas*) 

1 

0.14 

0.20 

0.14 

0.7 

2 

0.62 

0.36 

0.14 

0.8 

3 

1.77 

1.04 

0.24 

1.8 

4 

0.44 

0.42 

0.16 

1.0 

5 

1.30 

0.91 

0.18 

1.5 

6 

0.48 

0.31 

0.14 

0.9 

7 

1.55 

1.03 

0.21 

1.6 

8 

0.51 

0.28 

0.15 

0.8 

9 

0.97 

0.74 

0.14 

1.2 

10 

0.67 

0.40 

0.14 

0.9 

Average 

0.84 

0.57 

0.16 

1.1 

* mas  is  the  unit  abbreviation  for  0".001. 


deviation  in  the  parallax  of  0".0005,  an  improvement  by  a factor  of  two  over  the  value 
given  in  table  4-2. 

Tables  4-3  and  4-4  present  similar  results  for  the  adjustments  incorporating  the 
"observations"  from  the  association  models  constructed  with  uniform  distributions.  Once 
again,  the  gross  parameters  of  the  associations  are  recovered  nicely,  despite  the  presence 
of  two  models  (numbers  5 and  6)  that  seem  to  be  statistical  "flukes".  These  models  exert 
a strong  influence  on  the  average  absolute  deviations,  causing  them  to  be  higher  than  their 
Gaussian  counterparts  above.  If  the  averages  are  formed  without  the  values  belonging  to 
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uniform  associations  5 and  6,  one  obtains  results  almost  identical  to  the  numbers  in  the  last 
line  of  table  4-3.  Judging  by  the  numbers  in  tables  4-1,  4-2,  4-3,  and  4-4,  it  would  appear 
that  the  adjustment  with  stochastic  constraints  is  not  terribly  sensitive  to  the  forms  of  the 
association  distance  and  velocity  distribution  functions,  provided  that  suitable  constraints 
are  chosen. 


Table  4-3.  Global  parameters  for  uniform  association  models 


Model  No. 

v„  (km  s-1) 

v„,  (km  s1) 

v«.  (km  s1) 

7t„  (arcsec.) 

1 

-1.91  ± 1.67 

44.70  ± 1.06 

6.30  ± 0.32 

0.0193  ± 0.0020 

2 

-0.16  ± 1.40 

43.94  ± 0.90 

6.79  ± 0.26 

0.0214  ± 0.0020 

3 

-2.30  ± 1.47 

44.84  ± 0.94 

6.36  ± 0.28 

0.0191  ±0.0018 

4 

-1.91  ± 1.67 

44.70  ± 1.06 

6.30  ± 0.32 

0.0193  ± 0.0020 

5 

-4.55  ± 1.82 

46.40  ± 1.16 

5.99  ± 0.35 

0.0172  ±0.0018 

6 

0.99  ± 1.31 

43.06  ±0.81 

6.74  ± 0.25 

0.0224  ± 0.0020 

7 

-2.38  ± 1.51 

45.25  ± 0.97 

6.22  ± 0.30 

0.0191  ±0.0018 

8 

-0.75  ± 1.51 

44.06  ± 0.96 

6.48  ± 0.28 

0.0209  ± 0.0020 

9 

-0.37  ± 1.43 

43.77  ±0.91 

6.56  ± 0.27 

0.0210  ± 0.0020 

10 

-0.50  ± 1.44 

44.01  ±0.92 

6.52  ± 0.28 

0.0209  ± 0.0020 

Average 

-1.38  ±1.55 

44.47  ± 0.93 

6.43  ± 0.24 

0.0201  ±0.0015 

Table  4-4. 

Deviations  of  stochastic  parameters  from  true  values  for  uniform  models. 

Model  No. 

1 Av„  1 (km  s-1) 

1 Avv  1 (km  s_1) 

1 Av7  1 (km  S'1) 

1 Arc  1 (mas) 

1 

0.50 

0.24 

0.22 

0.9 

2 

1.25 

0.57 

0.30 

1.3 

3 

0.89 

0.34 

0.16 

0.9 

4 

0.50 

0.24 

0.22 

0.9 

5 

3.14 

1.92 

0.50 

2.8 

6 

2.40 

1.47 

0.26 

2.4 

7 

0.97 

0.73 

0.28 

1.1 

8 

0.66 

0.46 

0.13 

1.0 

9 

1.04 

0.76 

0.16 

1.1 

10 

0.91 

0.51 

0.14 

1.1 

Average 

1.23 

0.72 

0.24 

1.4 

The  adjustments  involving  the  last  group  of  four  models  with  a different  cluster  center 
and  velocity  yielded  no  surprises,  demonstrating  that  the  adjustment  technique  is  not 
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sensitive  to  either  the  association's  location  or  the  direction  of  its  average  velocity. 
Therefore,  for  the  sake  of  brevity,  a summary  of  these  models  shall  not  be  presented. 

Comparison  with  Convergent  Point  Solutions 

In  chapter  2,  much  ado  was  made  of  the  fact  that  an  adjustment  by  means  of  least 
squares  with  probabilistic  constraints  should  yield  more  accurate  parameter  estimates  than 
other  techniques,  as  it  incorporates  all  kinematical  observations  and  relaxes  the  unrealistic 
assumption  of  a shared  velocity.  The  model  associations  with  the  Gaussian  distributions 
presented  a perfect  opportunity  to  test  this  statement,  and  so  some  of  the  ten  models  in  the 
first  group  were  incorporated  into  two  more  adjustments.  Parameters  for  models  1,  2,  4, 

6 and  8 were  first  estimated  by  the  stochastic  constraint  technique,  with  the  parallax 
constraint  tightened  to  a value  of  0".002  (a  5 parsec  dispersion).  Then  a solution  was 
obtained  for  each  of  the  models  by  means  of  a classical  convergent  point  technique,  similar 
to  that  used  by  Wayman  et  al.  (1965).  Finally,  the  estimates  obtained  by  each  of  the 
methods  were  compared  to  the  true  values  to  determine  their  accuracy. 

Table  4-5  is  a summary  of  the  stochastic  constraint  adjustments;  note  that  the  global 
velocity  parameters  are  not  significantly  different  from  the  corresponding  values  in  table 
4-1.  The  most  important  values  for  the  current  discussion  are  the  average  deviations  from 
the  true  parallaxes,  Ait,  and  the  average  absolute  deviation,  IArc|.  A glance  at  table  4-5  and 
table  4-6  is  sufficient  to  show  that  the  stochastic  constraint  adjustments  yield  more 
accurate  parallaxes  than  those  estimated  by  the  convergent  point  method  (43%  in  the  case 
of  |A7t|  and  100%  with  regard  to  Ati).  A comparison  of  the  individual  stellar  velocities  is 
somewhat  meaningless,  as  the  convergent  point  solution  yields  only  one  velocity,  shared 
by  all  stars  in  the  association. 
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Table  4-5.  Gaussian  associations  with  parallax  constraints  of  2 mas. 


Model 

No. 

vxc  (km  s-b 

vvc  (km  s-b 

v7C  (km  s-b 

7ic  (mas) 

An 

(mas) 

1 Aji  1 
(mas) 

1 

-1.3±1.2 

44.610.8 

6.5  1 0.2 

1991 1.6 

.03 

.4 

2 

-0.5  ± 1.2 

44.0 1 0.8 

6.5  1 0.3 

208  1 1.7 

.8 

.9 

4 

-1.81 1.2 

44.9  1 0.8 

6.4  1 0.2 

1951 1.5 

-.6 

.7 

6 

-0.81 1.1 

44.2 1 0.7 

6.6  1 0.2 

204 1 1.6 

.4 

.7 

8 

-0.91 1.1 

44.3  1 0.7 

6.4 1 0.2 

204 1 1.6 

.3 

.6 

Average 

-1.1 

44.4 

6.5 

202 

.21.5 

.71.2 

Table  4-6. 
Model  No. 

Classical  convergent  point  solutions, 
a,  5 CP  Atc  (mas)  1 An  1 (mas) 

1 

91°.7,  8°.2 

.1 

.8 

2 

90°.7,  8°.4 

.9 

1.2 

4 

91°.8,  8°. 2 

-.08 

.8 

6 

91°.5,  8°.3 

.2 

.9 

8 

90°.8,  8°.4 

.8 

1.1 

Average 

91°.3,  8°.3 

.41.4 

1.010.2 

More  evidence  of  the  superiority  of  the  stochastic  constraint  technique  is  provided  by 
figure  4-5,  which  shows  the  distributions  of  the  true  parallaxes,  the  parallaxes  obtained  by 
the  stochastic  adjustment,  and  the  convergent  point  parallaxes  for  one  of  the  model 
associations.  The  stochastic  parallax  distribution  (figure  4-5b)  is  a fair  approximation  to 
that  of  the  true  parallaxes  (figure  4-5a),  but  the  convergent  point  parallax  distribution 
(figure  4-5c)  is  much  wider  than  that  of  the  true  values.  The  stochastic  adjustment  has 
preserved  the  shape  of  the  association,  whereas  the  convergent  point  parallaxes,  as 
expected,  "smear  out"  the  stars,  distorting  the  association  along  the  line  of  sight.  Tables 
4-5  and  4-6,  combined  with  figure  4-5,  should  provide  sufficient  reason  for  anyone 
concerned  with  accuracy  to  abandon  the  determination  of  parallaxes  by  means  of  the 
classical  convergent  point  approach. 

However,  it  would  be  misleading  to  state  that  parameters  obtained  by  an  adjustment 
with  stochastic  constraints  will  always  be  accurate  estimators  of  the  true  values.  As  with 


Number  of  stars  Number  of  stars  Number  of  stars 
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Figure  4-5.  Distribution  of  model  parallaxes, 
a)  "True"  parallaxes  of  model  stars;  b)  Parallaxes  of  model  stars  obtained 
by  stochastic  constraint  technique;  c)  Parallaxes  of  model  stars 
obtained  by  classical  convergent  point  solution. 


67 


any  least-squares  procedure,  the  technique  assumes  that  the  observations  are  bias-free;  use 
of  biased  observations  will  propagate  the  bias  into  the  parameter  estimates.  Figure  4-6 
illustrates  the  effect  of  such  a bias  upon  the  solution.  The  parallax  estimates  obtained  by 
the  adjustment  (figure  4-6c)  again  reproduce  the  form  of  the  true  distribution  (figure 
4-6a),  but  are  offset  by  approximately  -0".001.  As  all  the  "observed"  parallaxes  of  the 
model  stars  enter  into  the  adjustment  with  equal  weight  (being  from  the  same  distribution), 
it  is  thought  that  the  near- zero  parallaxes  are  responsible  for  the  shift  of  the  association's 
center  from  50  parsecs  to  about  53  parsecs.  The  fault  lies  not  with  the  specific  technique 
of  probabilistic  constraints;  rather,  it  lies  with  the  inability  of  least  squares  to  handle  biased 
data.3  The  lesson  here  is  clear;  use  of  biased  observations  in  a least-squares  adjustment 
should  be  avoided  if  at  all  possible,  unless  the  bias  can  be  removed  by  some  means.  This 
will  have  an  important  bearing  on  the  choice  of  the  data  to  be  used  in  the  determination  of 
the  Hyades  distance. 


Application  of  the  Technique  to  the  Hvades 
Validity  of  Assumptions  and  Determination  of  Constraints 

In  the  adjustment  procedure  we  will  stochastically  constrain  the  velocity  components 
and  the  parallaxes  of  the  Hyades  stars  to  obey  Gaussian  distributions.  It  has  been  shown 
previously  in  this  chapter  that  the  parameter  estimates  obtained  by  the  least  squares  with 
probabilistic  constraints  are  not  terribly  sensitive  to  the  forms  of  these  distributions;  still,  it 
would  be  reassuring  if  some  idea  of  the  actual  distributions  could  be  obtained  before 
performing  the  adjustment.  Fortunately,  there  exist  excellent  dynamical  studies  of  the 


3 This  bias  originates  in  the  fact  that  randomly  distributed  small  parallaxes  introduce  a bias  into  the 

dn 

distances  because  dr  = -777;  other  things  being  equal,  a change  in  a small  parallax  produces  a larger 
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change  in  the  corresponding  distance  than  the  same  change  in  a large  parallax. 
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Hyades  by  van  Bueren  (1952),  Oort  and  Pels  (1975),  Gunn  et  al.  (1988)  and  Reid  (1992) 
that  shed  considerable  light  on  the  nature  of  the  distance  and  velocity  distributions. 

Oort  and  Pels  (1975)  found  that  the  surface  density  of  the  Hyades  stars  drops  off 
exponentially  with  distance  from  the  center  of  the  cluster.  Figure  4-7  is  a comparison  of 
their  surface  density  distribution  with  that  of  Reid  (1992),  who,  along  with  Gunn  et  al. 
(1988),  found  that  the  surface  density  distribution  of  the  Hyades  stars  could  be 
characterized  by  a Plummer  model.  In  this  simple  polytropic  approximation,  the  surface 
density  of  the  stars  is  given  by 


L(r)  = A 


( r2  V2 

1+- 7 

v A J 


(4-11) 


where  r is  the  angular  distance  in  degrees  from  the  center  of  the  cluster  and  N and  A are 
constants.  Note  that  in  the  case  of  the  Hyades,  the  Plummer  model  and  the  exponential 
model  of  Oort  and  Pels  give  approximately  the  same  results,  which  is  not  too  surprising  as 
the  cluster  stars  show  a concentration  towards  the  center.  It  should  be  pointed  out  that 
the  distributions  shown  in  figure  4-7  are  for  the  brighter  members  of  the  Hyades,  as  the 
Hyades  dwarfs  are  more  dispersed  and  have  a lesser  affinity  for  the  center.  As  a result, 
Oort  and  Pels  give  two  exponential  spatial  density  distributions,  one  for  the  dwarfs  and  the 
other  for  the  brighter  stars  (mainly  of  spectral  type  A).  Figure  4-8  shows  their  observed 
number  distribution  function  for  the  brighter  Hyads.  It  should  now  be  obvious  that  the 
distance  distribution  of  the  Hyades  stars  is  anything  but  uniform,  and  that  a Gaussian  with 
a suitably  chosen  width  can  approximate  this  distribution. 

An  approximate  width  of  the  distance  distribution  can  be  determined  from  figure  4-8, 
which  would  seem  to  indicate  a value  of  about  6 parsecs.  This  is  borne  out  by  an 
examination  of  figure  4-9,  which  is  a plot  of  the  equatorial  positions  of  the  Hyades  stars 
with  visual  magnitudes  brighter  than  9.5.  Even  though  the  cluster  is  of  large  extent 
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Figure  4-6.  Model  solution  showing  bias  in  stochastic  parallaxes, 
a)  True  parallaxes  of  model  stars;  b)  "Observed"  parallaxes;  c)  Stochastic  parallaxes. 
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gure  4-7. 


Normalized  surface  density  distribution  for  Hyades  stars. 


Distance  from  center  (pc) 

Figure  4-8.  Observed  distance  distribution  for  the  brighter  Hyads. 
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Figure  4-9.  Equatorial  positions  of  Hyades  stars  with  mv  > 9 (orientation  B 1950.0). 
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(>15°),  there  is  an  obvious  concentration  of  stars  in  the  area  defined  by  63°  < a < 73°, 

13°  < 5 < 23°,  or  for  a radius  of  about  5°  from  the  center  of  the  cluster.  Assuming  a 
Hyades  distance  of  48  parsecs,  5°  corresponds  to  a distance  of  4.2  parsecs  from  the 
center,  in  rough  agreement  with  the  value  determined  from  figure  4-8.  Given  these  values, 
and  also  recalling  that  the  faint  members  of  the  cluster  show  less  of  a concentration 
towards  the  center,  6 parsecs  as  our  distance  dispersion  appears  to  be  a reasonable  choice, 
resulting  in  a parallax  constraint  of  approximately  0".003.  It  should  be  pointed  out  that 
even  if  we  adopt  the  Hyades  tidal  limit  of  approximately  10  parsecs  (see  Oort  and  Pels, 
1975  or  Reid,  1992)  as  the  dispersion  width,  the  parallax  constraint  would  increase  by 
0".002,  to  a value  of0".005. 

Obtaining  the  dispersion  of  the  stars'  velocity  components  is  much  easier,  thanks  to  the 
radial  velocity  studies.  In  their  study  of  the  Hyades,  Oort  and  Pels  (1975)  assumed  that 
the  cluster  is  in  dynamical  equilibrium,  with  the  internal  motions  of  the  member  stars 
obeying  a (truncated)  Maxwellian  velocity  distribution.  Using  radial  velocities  from  the 
Wilson  catalogue,  they  estimated  a velocity  dispersion  of  approximately  0.25  km  s4  in 
each  cartesian  direction.  Gunn  et  al.  (1988)  assumed  a velocity  distribution  dictated  by 
the  adoption  of  the  Plummer  model  and  derived  a dispersion  value  of  0.23  km  s4,  which 
they  estimate  cannot  be  in  error  by  more  than  20%.  It  is  hard  to  argue  with  such  good 
agreement;  therefore,  we  shall  adopt  0.25  km  s4  as  the  stochastic  constraint  for  the 
velocities. 

Choice  of  Data 

The  model  associations,  while  demonstrating  the  superiority  of  the  technique  of  least 
squares  with  probabilistic  constraints  over  the  conventional  techniques,  also  illustrate  the 
dangers  of  using  biased  observations  in  the  adjustment.  Biases  are  easy  to  discover  in  the 
virtual  reality  of  computer  models,  but  nature  restricts  us  to  obtaining  only  estimates  of 
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the  quantities  we  wish  to  observe;  the  "true"  values  are  forever  out  of  reach,  cloaked  by 
the  unavoidable  errors  in  the  observing  process  and  the  equally-unavoidable  imperfections 
of  the  mathematical  reduction  models.  It  is  hoped,  and  all  too  often  assumed,  that  these 
errors  are  random  in  nature,  and  that  as  the  precision  of  the  observations  increases,  so  to 
does  their  accuracy.  Unfortunately,  the  history  of  astrometry  is  replete  with  examples  of 
estimates  based  on  inaccurate,  but  precise,  observations.  One  such  example  involves 
Hyad  proper  motion  determinations  made  by  various  investigators;  an  examination  of  each 
investigators'  statements  concerning  the  accuracies  of  their  motions  relative  to  those  of 
their  colleagues  evokes  the  feeling  that  one  has  just  chanced  upon  a raging  battle.  Indeed, 
the  fury  of  the  criticisms  (most  of  which  have  at  least  a degree  of  validity)  tends  to  cast 
doubt  as  whether  there  exist  any  proper  motions  reasonably  free  of  systematic  effects.4 

Nonetheless,  proper  motions  are  a necessary  element  in  determining  the  Hyades 
distance.  They  must  be  incorporated  into  the  adjustment,  requiring  the  development  of 
some  sort  of  selection  procedure.  One  possible  approach  is  to  combine  all  the  available 
Hyades  proper  motions  into  one  huge  data  set  and  hope  that  the  systematic  errors  will 
"cancel  each  other  out".  However,  these  proper  motions  are  a bewildering  mess  of 
differing  measurement  techniques  and  coordinate  systems,  thus  rendering  such  an 
adjustment  impracticable,  if  not  impossible.  Another  approach  would  be  to  adopt  the 
concept  of  "newer  is  better"  and  consider  only  the  most  modem  proper  motions  in  the 
adjustment,  justified,  perhaps  to  some  extent  by  the  consideration  that  proper  motions 
with  respect  to  remote  galaxies  are  probably  less  affected  by  systematic  errors  than  those 
based  on  transit  circle  observations  and  other  differential  photographic  techniques.  This 
approach  is  followed  by  most  investigators,  and  while  the  improvements  in  instruments 


4 Eichhom  (1993,  private  communication)  relates  that  at  the  First  Cosmic  Distance  Conference,  held  in 
1955  at  the  University  of  Virginia  in  Charlottesville,  W.  J.  Luyten  began  an  invited  talk  saying:  “I  am 
supposed  to  talk  about  proper  motions.  Unfortunately,  there  exist  only  improper  motions. . .”. 
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and  advances  in  data  acquisition  give  it  some  credence,  it  does  a disservice  to  the  older 
proper  motion  estimates.  Even  more  importantly,  the  information  provided  by  these  older 
motions  is  ignored.  It  would  seem  that  there  is  no  good  approach,  the  choice  being 
between  loss  of  information  by  ignoring  past  data  or  reliance  on  the  chance  cancellation  of 
the  systematic  errors  of  the  various  proper  motions  in  a global  adjustment. 

However,  it  is  possible  to  compare  the  proper  motions  of  stars  common  to  multiple 
investigations;  hopefully,  such  a comparison  will  provide  some  insight  as  to  which  proper 
motions  suffer  from  bias  the  least.  It  has  long  been  known  that  the  transit  circle  proper 
motions  used  in  the  older  Hyades  distance  determinations  suffer  from  regional  systematic 
errors  that  limit  their  usefulness  (Hanson,  1975;  Hanson,  1980;  Reid,  1992).  The  sole 
exception  seems  to  be  the  FK  5 proper  motions  determined  by  Schwan  (1992),  which 
yield  a Hyades  distance  compatible  with  that  obtained  by  the  photographic  proper  motions 
(Reid,  1992).  This  would  indeed  be  a fortunate  circumstance,  as  the  radial  velocities  of 
140  of  his  145  Hyades  members  have  been  determined  as  well  as  the  absolute  parallaxes  of 
37  of  these  stars  (explained  by  the  fact  that  the  Schwan  stars  have  visual  magnitudes 
brighter  than  9.5). 

Systematic  errors  plague  the  galaxy-referred  photographic  proper  motions  as  well. 
Hanson's  (1975)  proper  motions  were  touted  as  being  free  of  systematic  effects  due  to  the 
use  of  external  galaxies  as  reference  objects.  Two  years  later,  McAlister  (1977)  showed 
that  these  proper  motions  suffered  from  a magnitude  dependency,  and  a comparison  with 
Reid  (1992)  also  shows  a significant  discrepancy  between  his  estimates  and  those  of 
Hanson.  However,  the  Schmidt  plate  proper  motions  of  Luyten  et  al.  (1981)  and  Reid 
(1992)  agree  nicely  and  have  not  been  shown  to  suffer  from  major  systematic  effects  as  of 
this  date,  making  them  good  candidates  for  incorporation  into  an  adjustment. 
Unfortunately,  the  quoted  precision  of  these  motions  is  much  less  than  that  of  Hanson’s 
(1975)  or  even  that  of  Schwan’s  (1991);  it  would  seem  that  accuracy  has  been  improved 
at  the  cost  of  precision. 
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The  use  of  comparisons  seems  to  reduce  the  set  of  usable  proper  motion  observations 
to  those  of  Schwan  (1991),  Luyten  et  al.  (1982)  and  Reid  (1992).  The  older  proper 
motions  of  Smart  (1939),  Van  Bueren  (1952),  Corbin  et  al.  (1975)  and  Oort  and  Pels 
(1975)  are  not  good  candidates  for  use  in  adjustments,  not  because  they  are  older  and 
more  imprecise,  but  because  they  suffer  from  serious  systematic  effects  that  would  bias  the 
parameter  estimates.  Hanson's  (1975)  proper  motions  could  possibly  be  corrected  for  the 
magnitude  effect,  but  the  rewards  would  probably  be  small,  as  there  is  a large  degree  of 
overlap  between  his  data  and  that  of  Reid  (1992). 

There  are  not  many  trigonometric  parallaxes  available,  and  those  that  exist  for  Hyades 
members  are  incorporated  into  the  adjustment,  provided  they  were  reduced  to  an  absolute 
system  (van  Altena,  1974a).  To  a lesser  extent  this  is  also  true  of  the  radial  velocities, 
where  the  observations  come  from  three  sources:  Wilson’s  catalogue  (standard  error  of  ~5 
km  s'1),  the  radial  velocities  of  Detweiler  et  al  (1984)  (standard  error  of  -0.25  km  s1),  and 
the  radial  velocities  of  Griffin  et  al.  (1988)  (standard  error  of  -0.25  km  s1).  There  may  be 
a zero  point  error  (bias)  of  about  1 km  s4  in  the  high  precision  estimates  due  to  the 
difficulties  involved  in  the  calibration  procedure.  However,  this  is  not  certain,  making 
these  velocities  prefered  over  those  from  Wilson’s  catalogue  due  to  their  much  higher 
precision. 

Table  4-7  lists  the  observations  for  the  140  stars  chosen  for  the  adjustment.  The  first 
column  lists  the  Schwan  (1992)  number,  and  the  next  two  columns  give  the  star's  B1950.0 
position  (in  units  of  decimal  degrees).  Columns  4 through  7 list  the  proper  motion 
components  and  their  associated  standard  errors  (units  are  104  arcseconds).  The  radial 
velocity  data  are  contained  in  columns  8 and  9;  a G in  column  9 indicates  that  the  velocity 
is  from  Gunn's  (1988)  data,  whereas  a W indicates  that  value  is  from  Wilson’s  catalogue. 
No  errors  are  given  for  the  radial  velocities;  these  are  assumed  to  be  0.25  km  s4  for 
Gunn’s  estimates  and  5 km  s4  for  Wilson’s  values.  Columns  10  and  1 1 list  the  visual 
magnitudes  and  the  B-V  indices  for  the  stars.  Finally,  the  observed  parallax  (if  available) 
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and  its  standard  error  are  given  for  each  star  in  columns  12  and  13  (parallax  units  are  ", 
whereas  the  error  is  given  in  units  of  mas). 

The  Adjustment. 

Once  the  constraints  had  been  established  and  the  observations  selected,  an  adjustment 
was  performed  using  a FORTRAN  program  similar  to  that  given  in  appendix  A.  The 
starting  values  of  the  parameters  were  set  at 

x0  =-5 

y0  =45 

z0  = 5.5 
7t0  =0"  .0225 

Convergence  was  achieved  after  ten  iterations,  and  the  resulting  global  parameters  are 
presented  in  table  4-8.  Note  that  the  parallax  estimate  places  the  cluster  center  at 
49.8  ±1.1  parsecs,  a value  in  good  agreement  with  recent  determinations  (Schwan,  1991 
and  Reid,  1992).  The  residuals  are  plotted  in  figures  4-11  through  4-15  (they  are  listed  in 
appendix  B).  Note  that  there  are  no  obvious  systematic  errors  in  any  of  the  residuals  with 
respect  to  position  or  the  stars'  visual  magnitude. 

The  stochastic  parameters,  the  "gems  in  the  crown"  of  this  adjustment  procedure,  are 
listed  in  table  4-9.  The  velocity  components  have  standard  deviations  of  0.22,  0.31  and 
0.40  km  s1  in  the  equatorial  x,  y,  and  z directions,  respectively.  The  stochastic  parallaxes 
have  a standard  deviation  of  0".0019  (4.7  parsecs)  about  their  mean  of  0".0201,  with 
extrema  of  0".0301  (33.2  parsecs)  and  0".0151  (66.2  parsecs).  These  parallaxes  are 
depicted  in  figure  4-10. 
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23 
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26 

27 

28 
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Table  4-7.  Observations  chosen  for  incorporation  into  the  adjustment. 


a 

8 

ds  es 

rv  I 

m„  B-V 

7t 

64.2375 

15.5167 

1198  4 

-245  4 

39.30  G 

3.64  0.99 

0.0263 

100 

65.0125 

17.4333 

1122  4 

-291  5 

39.70  G 

3.76  0.98 

0.0188 

90 

66.4250 

19.0667 

1135  3 

-372  4 

39.40  G 

3.54  1.02 

0.0155 

131 

75.0250 

21.5167 

706  4 

-414  5 

42.20  W 

4.64  0.15 

0.0185 

44 

43.8042 

20.4667 

2488  6 

-298  6 

28.50  W 

5.81  0.41 

57.5750 

17.1833 

1485  7 

-298  7 

35.00  W 

5.97  0.34 

0.0310 

85 

67.7500 

14.7333 

1063  6 

-261  7 

37.50  W 

4.65  0.25 

0.0241 

85 

72.1125 

18.7500 

843  6 

-345  7 

38.50  W 

5.10  0.22 

0.0214 

80 

76.6542 

27.9667 

619  9 

-598  9 

45.00  G 

5.99  0.27 

76.6417 

9.7667 

645  7 

-71  8 

37.20  W 

5.42  0.24 

0.0099 

119 

62.1667 

5.4000 

1492  10 

55  11 

36.60  W 

5.72  0.36 

0.0320 

92 

63.8542 

21.4667 

1014  9 

-366  10 

35.40  W 

5.65  0.28 

65.8292 

22.7000 

1162  13 

-458  13 

35.10  W 

4.28  0.26 

0.0348 

170 

66.5083 

12.9333 

1074  19 

-156  22 

33.40  W 

5.03  0.22 

68.8417 

12.4167 

1057  16 

-128  19 

45.00  W 

4.27  0.12 

0.0216 

72 

69.1000 

7.7667 

948  19 

-29  22 

36.30  W 

5.39  0.26 

69.1042 

15.8167 

864  13 

-181  14 

36.90  W 

4.69  0.15 

70.8125 

11.6167 

691  19 

-52  23 

40.80  W 

5.37  0.20 

0.0043 

85 

59.2333 

10.1833 

1419  33 

-62  30 

39.60  W 

6.38  0.41 

63.2333 

15.2833 

1174  22 

-245  21 

36.90  W 

6.33  0.40 

64.2875 

13.9167 

1216  22 

-182  24 

42.00  W 

5.59  0.28 

64.4417 

14.9833 

1131  18 

-202  22 

36.20  W 

5.26  0.23 

64.5167 

13.7500 

1092  18 

-207  17 

37.00  W 

6.18  0.46 

64.8083 

13.9667 

1143  16 

-240  19 

41.20  W 

5.72  0.31 

64.9667 

16.6833 

1464  39 

-275  32 

33.00  W 

7.86  0.68 

65.1375 

16.6667 

1123  18 

-174  22 

35.00  W 

5.64  0.30 

65.3042 

17.3333 

1134  13 

-318  16 

37.50  W 

4.80  0.16 

0.0167 

136 

65.5083 

18.9333 

1125  16 

-345  16 

36.60  W 

5.98  0.37 

65.5958 

22.1833 

1135  19 

-407  19 

40.10  W 

4.22  0.14 

0.0229 

138 

65.6083 

22.0833 

1197  19 

-442  19 

32.00  W 

5.27  0.24 

65.6500 

17.8167 

1158  13 

-364  17 

40.20  G 

4.29  0.05 

0.0232 

28 

65.6917 

15.8333 

1050  31 

-253  35 

36.40  W 

6.46  0.49 

0.0370 

153 

66.2583 

21.5167 

1059  16 

-413  18 

36.20  W 

5.72  0.27 

66.4292 

15.8500 

1096  15 

-293  15 

40.20  G 

3.85  0.95 

0.0369 

155 

66.4500 

15.7667 

1117  16 

-255  18 

39.50  W 

3.40  0.18 

0.0272 

170 

66.9250 

16.0833 

1092  16 

-239  19 

37.50  W 

4.78  0.17 

0.0627 

101 

66.9500 

15.5833 

1072  12 

-247  18 

39.30  W 

5.48  0.26 

66.9500 

13.6167 

1135  21 

-193  25 

38.80  W 

5.40  0.26 

67.2500 

15.7500 

1075  24 

-295  23 

36.00  W 

6.02  0.34 

68.6708 

15.0500 

1095  19 

-292  18 

42.00  G 

7.54  0.61 

68.8208 

15.9333 

934  16 

-286  19 

38.40  W 

5.78  0.31 

68.9625 

23.0500 

1113  22 

-578  19 

39.90  G 

7.53  0.58 

70.4125 

11.0500 

984  21 

-119  26 

39.40  W 

5.40  0.25 

74.2167 

15.8500 

876  15 

-322  16 

44.10  W 

6.75  0.45 

52.4708 

23.5167 

1755  39 

-353  37 

32.70  G 

8.89  0.84 

53.6667 

21.1833 

1635  36 

-354  33 

33.60  G 

9.37  0.92 

62.8833 

22.3333 

1035  21 

-265  20 

38.40  W 

7.05  0.42 

62.8833 

23.4500 

1351  36 

-514  33 

38.20  G 

8.08  0.66 

63.1917 

20.7000 

1216  21 

-355  21 

38.50  G 

7.84  0.71 

Table  4-7  - continued 
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No. 

a 

8 

Mw.  e„ 

4s  es 

rv 

I 

nv  B-V 

n e, 

50 

63.3958 

21.7833 

1171  21 

-425  19 

38.50 

G 

9.15  0.82 

51 

64.6417 

20.9167 

1126  21 

-365  20 

37.10  W 

6.80  0.44 

52 

65.1042 

21.2667 

1081  39 

-405  38 

39.20 

G 

8.85  0.76 

53 

65.1208 

24.2833 

1021  24 

-455  22 

34.10  W 

7.18  0.45 

54 

65.3167 

21.6167 

1276  21 

-484  21 

38.90 

G 

7.13  0.52 

55 

65.8750 

21.3667 

1111  21 

-345  20 

38.80 

G 

7.38  0.54 

56 

66.6083 

26.5667 

991  27 

-516  23 

39.00 

G 

9.19  0.76 

57 

67.0792 

20.0333 

1066  39 

-345  38 

40.60 

G 

8.59  0.74 

58 

68.3667 

23.2333 

1141  25 

-495  24 

43.30  W 

6.05  0.39 

59 

69.2708 

23.2000 

721  40 

-327  37 

41.30 

G 

9.40  0.81 

60 

71.4333 

21.0167 

916  19 

-436  18 

41.80 

G 

9.07  0.84 

61 

71.5417 

24.7167 

993  43 

-555  39 

40.60 

G 

9.59  1.05 

62 

76.6542 

27.9667 

543  51 

-477  42 

45.00 

G 

8.59  0.69 

63 

58.0625 

16.8500 

1500  22 

-274  22 

36.80 

G 

8.99  0.90 

64 

60.8583 

15.5667 

1260  42 

-215  42 

38.00 

G 

7.85  0.59 

65 

61.2167 

15.0333 

1440  51 

-254  44 

36.40 

G 

6.02  0.40 

66 

61.9542 

18.3000 

1290  24 

-264  23 

33.00  W 

6.62  0.42 

67 

62.9000 

14.5000 

1215  42 

-235  38 

39.20 

G 

8.46  0.70 

68 

62.9167 

12.3167 

1185  24 

-55  22 

39.30 

G 

8.06  0.64 

69 

62.9542 

10.5833 

1200  24 

-155  23 

38.50 

G 

7.12  0.51 

70 

63.6958 

16.8333 

1125  22 

-265  21 

38.20 

G 

8.32  0.76 

71 

63.7833 

18.1333 

1170  19 

-295  20 

44.00  W 

7.53  0.68 

72 

63.8667 

15.9667 

1245  22 

-275  21 

39.60 

G 

9.59  0.98 

73 

64.0083 

19.7833 

1186  39 

-395  38 

39.60 

G 

8.63  0.74 

74 

64.0625 

17.4000 

1155  22 

-205  22 

38.80 

G 

8.43  0.73 

75 

64.2625 

16.4000 

1185  48 

-285  41 

39.20 

G 

6.89  0.56 

76 

64.3250 

19.1167 

1231  21 

-355  21 

38.80 

G 

7.47  0.57 

77 

64.3792 

18.6167 

1215  25 

-305  24 

42.00  W 

6.13  0.37 

0.0190 

78 

64.6583 

18.3000 

1170  22 

-235  21 

37.00  W 

6.81  0.44 

79 

64.6875 

14.3000 

1125  25 

-155  24 

39.20 

G 

6.61  0.41 

80 

64.9167 

11.1833 

1215  22 

-155  22 

40.80 

G 

9.81  1.18 

81 

64.9750 

14.9333 

1320  45 

-184  40 

42.70 

G 

6.99  0.56 

82 

65.1125 

19.5500 

1036  39 

-405  39 

39.90 

G 

9.40  0.91 

83 

65.1417 

15.6500 

1276  55 

-364  54 

39.30 

G 

10.50  1.23 

0.0304 

84 

65.1750 

14.5500 

1215  22 

-185  21 

40.30 

G 

9.11  0.86 

85 

65.3417 

14.6500 

1185  22 

-205  22 

37.50 

G 

7.62  0.60 

0.0197 

86 

65.3375 

16.2667 

1035  42 

-245  41 

39.80 

G 

8.24  0.59 

87 

65.3750 

16.9667 

1051  46 

-285  40 

43.20  W 

6.97  0.44 

88 

65.4000 

16.7667 

1260  42 

-254  41 

38.80 

G 

7.80  0.60 

89 

65.7375 

18.7500 

1201  36 

-315  35 

40.70 

G 

7.53  0.68 

90 

65.8125 

15.4167 

mi  24 

-285  22 

40.00 

G 

7.48  0.54 

91 

65.8750 

15.5000 

1170  48 

-235  47 

41.40  W 

4.49  0.25 

0.0221 

92 

65.8833 

16.7333 

1170  25 

-245  24 

39.40 

G 

8.03  0.65 

0.0317 

93 

65.9500 

16.6333 

1155  24 

-285  24 

40.20 

G 

8.10  0.67 

0.0200 

94 

66.0125 

10.7667 

1065  24 

-115  22 

41.10 

G 

9.46  1.03 

95 

66.1875 

15.4833 

1140  25 

-265  24 

39.30 

G 

7.42  0.54 

0.0321  1 

96 

66.2458 

11.6333 

1080  24 

-105  22 

40.20 

G 

7.51  0.55 

97 

66.2583 

18.9500 

856  27 

-316  26 

41.00 

G 

11.00  1.20 

98 

66.3125 

13.7500 

1005  48 

-126  45 

40.90 

G 

9.02  0.84 

99 

66.3875 

14.6333 

1095  28 

-195  26 

44.20  W 

5.89  0.32 

Table  4-7  — continued 
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No. 

a 

5 

e„ 

V-*  e« 

rv  I 

mv  B-V 

K 

100 

66.4208 

19.6333 

1051  37 

-345  38 

39.90  G 

8.63  0.75 

0.0254 

31 

101 

66.4792 

17.1833 

1171  22 

-315  21 

40.00  G 

7.85  0.60 

0.0458 

200 

102 

66.5333 

16.0500 

1020  43 

-225  38 

45.00  W 

6.59  0.53 

103 

66.6125 

17.4333 

1111  24 

-355  22 

39.90  G 

7.04  0.50 

104 

66.6542 

17.7500 

1111  28 

-355  29 

38.60  G 

6.91  0.45 

0.0309 

100 

105 

66.6583 

17.7833 

1126  21 

-325  20 

40.70  G 

8.96  0.83 

106 

66.7708 

16.5667 

1066  22 

-275  21 

41.00  G 

9.10  0.86 

107 

66.8417 

19.7333 

1021  19 

-325  20 

38.00  W 

7.10  0.47 

108 

66.9333 

15.6333 

975  24 

-236  22 

40.80  G 

8.94  0.85 

0.0214 

29 

109 

66.9792 

16.0500 

1066  27 

-275  26 

43.60  W 

6.51  0.43 

110 

67.0458 

10.6500 

930  25 

-116  23 

41.40  G 

7.04  0.47 

0.0171 

24 

111 

67.4917 

15.9000 

931  42 

-346  41 

43.00  W 

8.94  0.88 

0.0465 

99 

112 

67.5333 

15.7167 

946  22 

-266  21 

41.60  G 

8.65  0.74 

113 

67.6875 

16.6667 

841  42 

-296  40 

41.30  G 

9.41  0.89 

0.0326 

48 

114 

67.7417 

13.1500 

975  21 

-116  21 

35.00  W 

6.62  0.43 

0.0345 

46 

115 

67.7833 

15.0667 

991  19 

-296  17 

37.60  W 

8.49  0.84 

116 

67.9208 

15.7167 

811  24 

-236  22 

41.10  G 

9.67  0.92 

117 

67.9333 

15.4000 

961  22 

-306  22 

40.80  G 

7.89  0.64 

118 

68.3083 

15.5833 

856  42 

-326  41 

41.50  G 

9.38  0.86 

0.0184 

33 

119 

68.4542 

15.7667 

901  27 

-306  25 

33.70  W 

6.65  0.44 

120 

69.0333 

14.0000 

990  22 

-235  21 

41.90  G 

7.96  0.66 

121 

69.2625 

12.6333 

1020  24 

-155  23 

42.50  G 

10.04  1.08 

122 

69.3875 

16.4167 

916  49 

-336  46 

40.90  G 

9.47  1.10 

123 

70.9125 

15.3833 

886  22 

-326  21 

42.40  G 

8.32  0.67 

124 

71.1792 

18.1667 

916  22 

-406  22 

42.00  W 

8.53  0.72 

125 

71.5333 

18.5500 

961  21 

-426  20 

41.60  G 

8.99  0.83 

126 

71.6667 

15.8000 

886  22 

-286  21 

42.00  G 

7.74  0.58 

127 

71.8750 

17.1167 

1006  24 

-365  24 

41.70  W 

7.11  0.56 

128 

71.9292 

15.0000 

841  22 

-226  21 

42.30  G 

7.59  0.73 

129 

71.9833 

16.1333 

826  27 

-336  25 

42.70  G 

7.29  0.50 

130 

72.1083 

10.9833 

825  48 

-136  42 

42.40  G 

6.77  0.53 

131 

72.1333 

15.3500 

646  24 

-207  23 

42.90  G 

7.89  0.53 

132 

72.2542 

13.5667 

795  48 

-206  44 

38.80  W 

6.27  0.50 

133 

73.0083 

19.4000 

721  42 

-337  39 

35.00  W 

6.37  0.29 

134 

73.7500 

13.9333 

660  21 

-187  19 

43.70  G 

8.89  0.74 

0.0247 

29 

135 

75.0750 

13.6667 

645  48 

-137  44 

43.50  G 

9.00  0.76 

136 

48.6875 

7.4667 

1725  22 

-54  21 

32.40  G 

7.39  0.57 

0.0267 

37 

137 

65.5167 

4.5833 

1020  40 

55  38 

42.20  G 

8.93  0.76 

138 

67.3542 

5.3000 

1140  25 

115  24 

40.30  W 

6.38  0.41 

139 

71.0083 

8.9333 

1080  25 

-45  23 

42.50  G 

7.26  0.56 

140 

75.7458 

6.4000 

675  48 

-37  47 

42.90  G 

9.92  0.95 

Table  4-8.  Hyades  cluster  parameter  estimates. 

xc  (km  s'1)  yc  (km  s'1)  zg  (km  s'1)  tc„  (mas) 

-7.5010.41  47.11  10.23 5.2310.13 20.1 10.4 


k (arcsec.) 
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Star  number 


Figure  4-10.  Stochastic  parallax  estimates. 


O'yr1)  v ("  yr1)  v ("yr1) 
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mv 

Figure  4-11.  Residuals  in  |ia. 

a)  Residuals  in  |ia;  b)  Residuals  in  |ia  vs.  a;  c)  Residual  in  \ xa  vs.  ir^. 


(”  yr1)  V 6 ("  yr1)  v 5 ("  yr1) 
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Figure  4-12.  Residuals  in  (i8. 

a)  Residuals  in  |i6;  b)  Residuals  in  (i6  vs.  8;  c)  Residuals  in  (i6  vs.  rr^. 
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Star  number 

Figure  4-13.  Residuals  in  the  radial  velocity  and  parallax, 
a)  Residuals  in  the  radial  velocity;  b)  Residuals  in  7C. 


Residual  in  vr  (km  s4)  Residual  in  vr  (km  s'1)  Residual  in  vr  (km  s4) 
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mv 

Figure  4-14.  Residuals  in  the  radial  velocity, 
a)  Residuals  in  vr  vs.  a;  b)  Residuals  in  vr  vs.  5;  c)  Residuals  in  vr  vs.  mv. 


(arcsec.)  (arcsec.)  (arcsec.) 
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Declination  (degrees) 


Figure  4-15.  Residuals  in  the  parallax, 
a)  Residuals  in  7t  vs.  a;  b)  Residuals  in  7t  vs.  8;  c)  Residuals  in  n vs.  nv 
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Table  4-9.  Stochastic  velocities  and  parallaxes  for  140  Hyades  members. 


No. 

x (km  s1) 

y (km  s1) 

z (km  s1) 

7t  (arcsec.) 

1 

-7.1070 

47.0220 

5.2040 

0.0204 

2 

-7.0690 

47.1170 

5.3160 

0.0201 

3 

-7.1580 

46.8440 

5.0770 

0.0203 

4 

-6.9650 

47.0570 

5.3010 

0.0180 

5 

-7.2800 

47.3410 

5.2840 

0.0300 

6 

-6.8930 

47.2210 

4.6350 

0.0211 

7 

-6.9040 

47.1710 

4.6820 

0.0203 

8 

-7.0630 

47.0540 

5.2970 

0.0195 

9 

-6.2250 

48.5490 

5.7550 

0.0171 

10 

-7.0620 

46.9910 

5.5090 

0.0194 

11 

-7.1300 

47.1900 

4.6930 

0.0219 

12 

-6.8860 

47.0690 

5.0680 

0.0178 

13 

-7.2470 

47.0210 

5.5250 

0.0202 

14 

-7.1030 

47.0290 

5.3430 

0.0201 

15 

-7.1800 

47.0760 

5.7590 

0.0206 

16 

-7.0130 

47.0600 

5.0920 

0.0192 

17 

-7.2090 

46.8900 

6.0050 

0.0173 

18 

-7.0530 

46.9740 

5.7980 

0.0178 

19 

-7.1060 

47.1710 

5.2510 

0.0220 

20 

-6.9840 

47.1370 

4.9720 

0.0202 

21 

-7.1290 

47.1330 

5.4010 

0.0213 

22 

-7.1060 

47.0540 

5.3970 

0.0200 

23 

-6.9270 

47.1370 

4.8170 

0.0195 

24 

-6.9080 

47.2190 

4.6650 

0.0205 

25 

-7.4470 

47.0830 

5.6580 

0.0251 

26 

-7.4280 

46.8970 

6.2590 

0.0199 

27 

-6.9620 

47.1380 

4.9700 

0.0202 

28 

-7.0860 

47.0890 

5.2510 

0.0204 

29 

-7.2940 

47.0350 

5.6820 

0.0201 

30 

-7.2970 

47.0270 

5.4910 

0.0216 

31 

-6.8330 

47.3560 

4.7730 

0.0212 

32 

-7.0290 

47.0910 

5.1450 

0.0200 

33 

-7.0430 

47.0760 

5.2240 

0.0197 

34 

-6.9360 

47.0790 

4.7180 

0.0202 

35 

-7.0840 

47.0920 

5.3190 

0.0202 

36 

-7.2020 

47.0240 

5.6160 

0.0202 

37 

-7.0810 

47.0940 

5.2930 

0.0203 

38 

-7.1360 

47.1060 

5.3170 

0.0214 

39 

-6.9750 

47.1350 

4.8860 

0.0206 

40 

-6.9090 

47.4490 

4.9100 

0.0218 

41 

-6.8510 

47.1310 

4.8180 

0.0189 

42 

-6.9850 

47.0380 

4.7340 

0.0225 

43 

-7.1000 

47.0840 

5.3060 

0.0206 

44 

-7.0070 

47.1790 

4.9970 

0.0211 

45 

-7.2820 

47.1730 

5.6150 

0.0236 

46 

-7.1100 

47.1920 

5.2360 

0.0227 

47 

-7.5070 

46.8490 

6.3270 

0.0172 

48 

-7.0500 

47.3390 

5.2880 

0.0227 

49 

-7.2150 

47.1020 

5.6030 

0.0206 

Table  4-9  — continued 
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No. 

x (km  s_1) 

y (km  S'1) 

z (km  s1) 

7i  (arcsec.) 

50 

-6.9320 

47.2230 

5.1080 

0.0202 

51 

-7.1840 

47.0410 

5.5050 

0.0199 

52 

-6.9910 

47.1330 

5.1960 

0.0196 

53 

-6.9930 

47.0450 

5.2110 

0.0185 

54 

-7.1850 

47.0660 

5.1060 

0.0229 

55 

-7.5760 

46.6830 

5.8860 

0.0199 

56 

-7.0120 

47.1800 

5.4170 

0.0186 

57 

-7.0620 

47.2580 

5.4990 

0.0200 

58 

-7.3530 

47.1030 

5.5990 

0.0219 

59 

-6.8190 

47.3280 

5.6360 

0.0152 

60 

-7.0240 

47.1640 

5.2620 

0.0199 

61 

-7.2440 

46.9500 

5.2800 

0.0210 

62 

-6.6170 

48.3940 

6.2090 

0.0154 

63 

-7.0060 

47.3520 

5.2000 

0.0228 

64 

-7.0450 

47.1480 

5.2630 

0.0204 

65 

-7.5380 

46.4500 

4.9160 

0.0230 

66 

-7.3170 

47.0020 

5.7510 

0.0212 

67 

-6.9780 

47.2390 

5.0840 

0.0207 

68 

-7.1880 

47.0520 

6.1200 

0.0200 

69 

-7.1200 

46.8920 

4.4260 

0.0205 

70 

-7.2130 

46.7040 

5.0620 

0.0196 

71 

-7.1030 

47.1290 

5.4240 

0.0203 

72 

-7.0070 

47.2900 

5.1420 

0.0216 

73 

-6.8680 

47.4470 

5.2100 

0.0208 

74 

-7.4530 

46.7750 

6.0300 

0.0199 

75 

-7.1010 

47.0340 

5.1200 

0.0208 

76 

-7.2200 

46.9810 

5.2750 

0.0215 

77 

-7.2210 

47.0880 

5.6100 

0.0208 

78 

-7.4650 

46.9300 

6.2000 

0.0203 

79 

-7.2890 

46.7640 

5.5880 

0.0199 

80 

-6.9060 

47.4960 

4.8830 

0.0217 

81 

-6.6940 

48.3380 

5.9250 

0.0228 

82 

-6.7610 

47.3580 

4.9830 

0.0192 

83 

-7.1370 

46.9090 

4.8860 

0.0213 

84 

-7.1950 

47.1880 

5.7410 

0.0216 

85 

-7.6890 

45.9570 

5.1420 

0.0212 

86 

-7.0270 

47.0280 

5.2210 

0.0189 

87 

-6.9730 

47.1370 

5.1740 

0.0193 

88 

-7.4510 

46.6310 

5.3210 

0.0222 

89 

-7.0550 

47.4790 

5.6520 

0.0215 

90 

-6.9680 

47.0900 

4.7530 

0.0205 

91 

-7.1070 

47.1290 

5.3130 

0.0210 

92 

-7.3910 

46.7280 

5.6130 

0.0207 

93 

-7.0900 

47.1320 

5.2950 

0.0207 

94 

-6.9090 

47.3330 

5.0930 

0.0197 

95 

-7.2370 

46.6520 

4.9930 

0.0203 

96 

-7.1900 

46.8490 

5.4260 

0.0200 

97 

-6.5630 

47.5260 

5.1070 

0.0165 

98 

-6.9880 

47.1880 

5.4900 

0.0188 

99 

-7.1090 

47.1130 

5.4730 

0.0203 

Table  4-9  — continued 
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No. 

x (km  S'1) 

y (km  s1) 

z (km  s1) 

7t  (arcsec.) 

100 

-7.0710 

47.0620 

5.2920 

0.0196 

101 

-7.1410 

46.9630 

5.2490 

0.0203 

102 

-7.0310 

47.1130 

5.3750 

0.0191 

103 

-7.0070 

47.0040 

4.7500 

0.0210 

104 

-7.2950 

46.4010 

4.7470 

0.0204 

105 

-7.0380 

47.2870 

5.3360 

0.0210 

106 

-6.9670 

47.3180 

5.3340 

0.0201 

107 

-7.1840 

47.0110 

5.5840 

0.0191 

108 

-6.9450 

47.1380 

5.2080 

0.0186 

109 

-7.0120 

47.1560 

5.1560 

0.0202 

110 

-6.8430 

47.2430 

4.9810 

0.0179 

111 

-6.8380 

47.2030 

4.7370 

0.0196 

112 

-6.7310 

47.4340 

4.9810 

0.0185 

113 

-6.7540 

47.2900 

5.0290 

0.0176 

114 

-7.1700 

46.9460 

5.8210 

0.0190 

115 

-6.7480 

47.2130 

4.3360 

0.0196 

116 

-6.7640 

47.0890 

4.9460 

0.0163 

117 

-6.8060 

47.0680 

4.4330 

0.0191 

118 

-6.7350 

47.2530 

4.7910 

0.0179 

119 

-6.7690 

47.1150 

4.5880 

0.0184 

120 

-6.9280 

47.2430 

4.9730 

0.0200 

121 

-7.0160 

47.3630 

5.4950 

0.0206 

122 

-7.0300 

46.8240 

4.8930 

0.0193 

123 

-6.6960 

47.2840 

4.4180 

0.0195 

124 

-6.7710 

47.2220 

4.6480 

0.0202 

125 

-6.9420 

46.9920 

4.6610 

0.0212 

126 

-7.0350 

46.8900 

4.9940 

0.0196 

127 

-7.1040 

47.1650 

5.1060 

0.0221 

128 

-7.0750 

46.8780 

5.3030 

0.0187 

129 

-6.7300 

47.2370 

4.6580 

0.0190 

130 

-7.0670 

46.7770 

5.0890 

0.0189 

131 

-6.7140 

47.1040 

5.1790 

0.0151 

132 

-6.9480 

47.0540 

5.1660 

0.0184 

133 

-6.9230 

46.9910 

5.2390 

0.0173 

134 

-6.7590 

47.1670 

5.2040 

0.0162 

135 

-6.9760 

46.8260 

5.4060 

0.0170 

136 

-6.6600 

47.7220 

4.0020 

0.0221 

137 

-6.6300 

47.8980 

5.1610 

0.0187 

138 

-7.1170 

47.0820 

5.7310 

0.0217 

139 

-7.2610 

47.0100 

5.4510 

0.0228 

140 

-7.1050 

46.3540 

4.9520 

0.0181 
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Conclusion 

This  aim  of  this  work  was  never  to  revise  the  cosmic  distance  scale  by  seeking  to 
determine  a Hyades  distance  radically  different  from  that  found  by  other  investigations. 
Such  a purpose  would  have  only  led  to  disappointment,  as  the  distance  estimate  of  49.7  ± 
1.0  parsecs  obtained  by  the  least  squares  with  probabilistic  constraints  is  in  good 
agreement  with  the  values  found  by  other  astronomers,  such  as  Reid  (1992).  Rather,  this 
thesis  had  as  its  goal  the  first-time  application  of  a new  adjustment  technique  to  a 
particular  problem  in  kinematic  astronomy,  one  that  has  piqued  the  interest  of  astronomers 
since  the  beginning  of  this  century.  This  procedure  is  capable  of  handling  all  types  of 
kinematic  data  — proper  motions,  radial  velocities  and  parallaxes;  no  information  is 
ignored.  It  enables  the  kinematics  of  the  cluster  to  be  incorporated  into  the  adjustment 
process,  thereby  eliminating  the  absurd  requirement  of  a common  velocity  for  all  stars 
which  are  members  of  the  Hyades.  The  parallaxes  of  these  stars  can  also  be  constrained  to 
prevent  the  observation  errors  from  distorting  the  cluster  along  the  line  of  sight.  Finally, 
it  has  now  been  demonstrated  that  this  adjustment  procedure  yields  more  accurate  velocity 
and  distance  estimates  than  the  standard  convergent  point  procedures  still  used  in  open 
cluster  investigations. 

The  only  excuses  that  may  be  given  for  not  using  a stochastically-constrained 
adjustment  to  estimate  the  kinematical  properties  of  open  clusters  are  that  it  requires  a 
good  deal  of  computing  power  and  considerable  time  spent  in  the  writing  and  debugging 
of  computer  programs.  These  excuses  seem  flimsy  in  light  of  the  efforts  spent  designing 
and  building  modem  astronomical  instruments.  Shall  we  spend  years  in  the  development 
and  construction  of  new  instrumentation,  only  to  see  the  data  provided  by  these 
instruments  wasted  in  adjustments  that  do  not  yield  the  best  possible  parameter  values? 
One  would  hope  not,  as  certainly  there  remains  much  work  to  be  done  before  the 
distances  to  the  Hyades  stars  can  be  definitively  determined,  especially  with  regard  to  the 
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biases  in  the  proper  motions.  It  has  been  remarked  that  the  data  provided  by  the 
HIPPARCOS  astrometric  satellite  will  provide  definitive  distance  estimates  of  all  known 
Hyades  members,  which  may  well  prove  to  be  true.  Whatever  the  quality  of  the 
HIPPARCOS  data,  this  much  is  certain;  the  distances  and  velocities  of  the  Hyades  stars 
have  not,  and  never  will  be,  definitively  established  by  any  convergent  point  solution. 
New  instruments  demand  new  adjustment  techniques. 


APPENDIX  A 

FORTRAN  ADJUSTMENT  PROGRAM 


Note:  This  program  does  not  check  for  convergence,  but  simply  performs  a set  number  of 
iterations,  outputting  the  values  of  the  parameter  corrections  on  each  pass. 

C234567 

c This  program  does  a constrained  least  square  fit  to  an  open 

cluster 

c 

c dimension  and  declare  arrays 

c 

implicit  real*8  (a-h,o-z) 
real*8  ks2pyr 
c 

c observation  vectors 

c 

dimension  ra(140) ,dc(140) ,pmra(140) ,pmdc(140) ,rv(140) ,plx(140) 
dimension  epmra (140) , epmdc (140) ,erv(140) , eplx (140) , nobsv (140) 
dimension  vmag(140) ,blessv(140) ,bdnum(140) 
c 

c initial  observation  vectors 

c 

dimension  raO (140) ,dc0 (140) ,pmra0 (140) ,pmdc0 (140) , rvO (140) 
dimension  plx0(140) 
c 

c major  array  dimensions 

c 

dimension  ah (420, 4) ,xh(420, 457) , sh(457, 457) , rhoh (560,560) 
dimension  w ( 420, 420) , atwa (4,4), bh (420,560) 
c 

c major  vector  dimensions 

c 

dimension  par (4) ,beta (560) , res (457) ,parcor (4) , fh (420) ,phi (420) 
dimension  pstar(140) ,vxstar(140) ,vystar(140) ,vzstar(140) 
dimension  respra (140) , respdc (140) ,resrv(140) , resplx (140) 
dimension  epar(4) 
c 

c work  array  dimensions 

c 

dimension  workl (457, 420) , work2 (420, 420) , work3 (420, 420) 
dimension  work4 (420, 4) , work5 (560, 420) 
c 

c work  vector  dimensions 

c 

dimension  wvecl (420) ,wvec2 (420) ,wvec3 (4) ,wvec4 (457) ,wvec5 (560) 
c 

c establish  constants  and  conversions 

c 
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pi=3 . 14159654d0 
dg2rad=pi/l . 8d2 
ks2pyr=l. 0277121 66D-6 
sc2rad=dg2rad/3 . 6d3 
c 

c set  up  dimensions 
c 

n=4 

m=140 

m3=420 

icx=457 

m4=560 

c 

c parameter  guesses 

c 

par (1) =-5 . 0d0*ks2pyr 
par (2) =4 . 50dl*ks2pyr 
par (3) =5 . 5d0*ks2pyr 
par (4) =2 ,25d-2 
c 

c constraints 

c 

sxdot=0 . 25*ks2pyr 
sydot=sxdot 
szdot=sxdot 
splx=2 . Od-3 
c 

c get  data 

c 

call  datard (bdnum, ra, dc, pmra, epmra, pmdc, epmdc, rv, 
& erv, plx, eplx, nobsv, vmag, blessv, m) 

c 
c 

c debug 

c 

write  (*,1999)  (nobsv(i),  i=l,40) 

1999  format (40i2) 
c 

c initialize  observation  vector 

c 

do  10  i=l,m 

raO (i) =ra (i) 
dcO (i) =dc (i) 
pmraO (i) =pmra (i) 
pmdcO ( i ) =pmdc ( i ) 
rvO (i) =rv (i) 

10  plxO (i) =plx (i) 
c 

c initialize  all  arrays 

c 

call  clearr (m3,n,ah) 
call  clearr (m3, icx, xh) 
call  clearr (m3, m4,bh) 
call  clearr (m3, m3, w) 
call  clearr (icx, icx, sh) 
call  clearr (m4,m4, rhoh) 
call  clearr (icx, m3, workl) 
call  clearr (m3, m3, work2) 
call  clearr (m3, m3, wo rk3) 
call  clearr (m3, n,work4) 
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call  clearr (m4, m3, work5 ) 

initialize  vectors 

call  clearv (icx, res) 
call  clearv(m4,beta) 
call  clearv (m4, wvec5) 
call  clearv ( icx, wvec4) 
call  clearv(m3, fh) 
call  clearv (m3, phi) 
call  clearv (m3, wvecl) 
call  clearv (m3, wvec2) 
call  clearv (n,wvec3) 
call  clearv (n,parcor) 

compute  sigma  and  rho  covariance  matrices 

call  getsig (epmra, epmdc, erv, eplx, nobsv, m, icx, sh) 
call  getrho (sxdot, sydot, szdot, splx,m4, rhoh) 

start  loop 

do  25  i=l, 10 

compute  a 

call  geta (par, ra, dc,m3, n, ah) 
calculate  x 

call  getx (par, ra, dc,m3, icx, nobsv, xh) 
write  (*,*)  'a  and  x computed' 

calculate  x * sigma  * x transpose 

call  mxmtg(sh,xh,icx,icx,m3,workl) 
call  mxmg (xh, workl,m3, icx, m3, work2) 
write  (*,*)  'x  * sigma  * xt  computed' 

compute  b 

call  getb(ah, m3, n,m4, nobsv, bh) 

calculate  b * rho  * b transpose 

call  mxmtg(rhoh,bh,m4,m4,m3,work5) 
call  mxmg(bh,work5,m3,m4,m3,work3) 
write  (*,*)  'b  * rho  * bt  computed' 

add  the  two 

call  matadd(work2, work3,m3,m3,w) 

invert  the  "weight"  matrix,  w 

call  spsinv(w,m3) 
write  (*,*)  'W  inverted' 


compute  phi  vector 
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call  mxvg(xh, res, m3, icx, wvecl) 

call  getf (ra,dc,par, rv,pmra,pmdc,m3, fh) 

call  vsubg(fh,wvecl,m3,phi) 

compute  corrections  to  the  parameters 

call  mxvg(w,phi,m3,m3,wvec2) 
call  mtxvg(ah,wvec2,n,m3,wvec3) 

compute  (a  transpose  * w * a)  **  -1 

call  mxmg(w,ah,m3,m3,n,work4) 
call  mtxmg (ah, work4, n,m3, n, atwa) 
call  syminv (atwa, n) 

bring  it  all  together 

call  mxvg (atwa, wvec3, n, n, parcor) 
call  negate (n, parcor) 

compute  the  residuals 

call  mxvg(ah, parcor, m3, n,wvecl) 
call  vaddg(phi,wvecl,m3,wvec2) 
call  mxvg (w, wvec2 ,m3,m3, wvecl ) 
call  mtxvg(xh, wvecl, icx, m3, wvec4) 
call  mxvg (sh,wvec4, icx, icx, res) 
call  negate (icx, res) 

correct  the  parameters 

call  vaddg (par, parcor, n, par) 

correct  the  observations  and  separate  the  residuals 
knt=l 

do  20  kk=l,m 
kkl=knt+l 
kk2=knt+2 
kk3=knt+3 

resrv (kk) =res (knt) /ks2pyr 
respra (kk) =res (kkl) /sc2rad 
respdc (kk) =res (kk2) /sc2rad 
resplx (kk) =res (kk3) 

rv (kk)  =rv0  (kk)  +res  (knt) 
pmra (kk) =pmra0 (kk) +res (kkl) 
pmdc (kk) =pmdc0 (kk) +res (kk2) 
plx (kk) =plx0 (kk) +res (kk3) 
knt=knt+nobsv (kk) 

20  continue 

outputs  to  check  for  convergence 

write  (*,600)  (parcor (j),  j=l,n) 

600  format (4e20 . 8) 


25  continue 
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c convergence  has  been  achieved,  compute  the  parameters 
c 

vxc=par ( 1 ) /ks2pyr 
vyc=par (2 ) /ks2pyr 
vzc=par (3) /ks2pyr 
plxc=par (4) 
c 

c now  calculate  the  beta  vector 
c 

call  mtxvg(bh,wvecl,m4,m3,wvec5) 
call  mxvg (rhoh, wvec5,m4,m4,beta) 
call  negate (m4, beta) 
c 

c compute  individual  star  velocities  and  parallaxes 
c 

do  30  j=l,m4,4 
knt= ( j-1) / 4+1 

jl-j+1 

j2=j+2 
j3= j+3 

vxstar (knt) = (par (1) +beta ( j) ) /ks2pyr 
vystar (knt) = (par (2) +beta ( jl) ) /ks2pyr 
vzstar (knt) = (par (3) +beta ( j2) ) /ks2pyr 
30  pstar (knt) =par (4) +beta ( j3) 
c 

c compute  errors  in  the  parameters 

c 

do  40  j=l,n 

epar ( j) =dsqrt (atwa ( j,  j) ) 
if  (j.lt.4)  then 

epar ( j ) =epar ( j ) /ks2pyr 
endif 
40  continue 
c 

c display  the  parameters  and  their  errors 
c 

write  (*,605)  vxc, vyc, vzc, plxc 
write  (*,605)  (epar(j),  j=l,4) 

605  format (3f8.2,f7. 4) 
c 

c open  output  file 

c 

open  (4, file='hyadesp .out ' ) 
c 

c write  results  to  output  file 

c 

do  50  j=l,m 

50  write  (4,  400)  resrv ( j ),  respra ( j ),  respdc (j ),  resplx ( j ), vxstar ( j ) , 
& vystar (j) , vzstar (j) , pstar (j) 

400  format (f7 .3, 3f 9 . 5, 3f 8 . 3, f 9 . 5) 
c 

c finished 

c 

close  (4) 

stop 

end 

C234567 

subroutine  getx (param, ra, dc,m3,nx, nobsv, x) 
c 

c INPUTS: 
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c param  is  the  parameter  vector:  x dot,  y dot,  z dot,  plx, 
respectively 

c ra  is  the  star  right  ascension  vector  (radians) 

c dc  is  the  star  declination  vector 

c m3  is  the  number  of  condition  equations 

c nx  is  the  column  dimension  of  x = sum  of  the  elements  of  nobsv 

c nobsv  is  the  vector  containing  the  number  of  observations  for  each 
c star  (2,  3 or  4) 

c 

c OUTPUTS: 

c x is  the  Jacobian  matrix  of  the  partials  of  the  condition 
equations 

c with  respect  to  the  observations 

c 

c declare  and  dimension  variables 

c 

implicit  real*8  (a-h,o-z) 

dimension  param (*) , ra  (*)  , dc  (*)  , nobsv  (*) , x (m3,  nx) 
c 

c initialize  column  pointer 

c 

jpos=l 

c 

c start  loop 

c 

do  50  i=l,m3,3 
knt= (i-1) /3+1 
i2=i+l 
i3=i+2 
j2=jpos+l 
j3= jpos+3 
j 4= jpos+2 
x (i, jpos) =1 . OdO 
x(i2,  j2) =1 . OdO 
x(i3, j4) =1 . OdO 
if  (nobsv (knt)  .gt.  3)  then 

x (i2,  j3) =param(l) *dsin (ra (knt) ) -param (2) *dcos (ra (knt) ) 
x (i3,  j 3 ) =param(l) *dcos (ra (knt) ) *dsin (dc (knt) ) + 

& param(2) *dsin (ra (knt) ) *dsin (dc (knt) ) - 

& param(3) *dcos (dc (knt) ) 

end  if 
c 

c increment  column  pointer 

c 

jpos= jpos+nobsv (knt) 

50  continue 
c 

c finished 

c 

end 

c 

subroutine  geta (param, ra, dc,m3, n, a) 
c 

c INPUTS: 

c param  is  the  parameter  vector:  x dot,  y dot,  z dot,  plx, 
respectively 

c ra  is  the  star  right  ascension  vector  (radians) 

c dc  is  the  star  declination  vector 

c m3  is  the  number  of  condition  equations 

c n is  number  of  unconstrained  parameters 
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c 

C OUTPUTS : 

c a is  the  Jacobian  matrix  of  partials  of  the  condition  equations 

with 

c respect  to  the  parameters 

c 

c declare  and  dimension  variables 

c 

implicit  real*8  (a-h,o-z) 
dimension  param(*) , ra (*) , dc ( *) , a (m3, n) 
c 

c start  loop 

c 

do  10  i=l,m3,3 
knt= (i-1) /3+1 
i2=i+l 
i3=i+2 

a (i,  1)  =-dcos  (ra  (knt) ) *dcos  (dc  (knt)  ) 
a (i,  2)  =-dsin  (ra  (knt) ) *dcos  (dc  (knt)  ) 
a (i,  3) =-dsin (dc (knt) ) 
a (i2, 1) =param(4) *dsin (ra (knt) ) 
a (i2, 2) =-param(4) *dcos (ra (knt) ) 

a (i2, 4) =param(l) *dsin (ra (knt) ) -param(2) *dcos (ra (knt) ) 
a (i3, 1) =param(4) *dcos (ra (knt) ) *dsin (dc (knt) ) 
a (i3, 2) =param( 4) *dsin (ra (knt) ) *dsin (dc (knt) ) 
a (i3, 3) =-param(4) *dcos (dc (knt) ) 
a (i3,  4)  =param(l)  *dcos  (ra  (knt)  ) *dsin  (dc  (knt)  ) + 

& param(2) *dsin (ra (knt) ) *dsin (dc (knt) ) - 

& param(3) *dcos (dc (knt) ) 

10  continue 
c 

c finished 

c 

return 

end 

c 

subroutine  getf  (ra,  dc,  param,  rv,pmra,pmdc,m3,  f) 
c 

c INPUTS : 

c param  is  the  parameter  vector:  x dot,  y dot,  z dot,  plx, 
respectively 

c ra  is  the  star  right  ascension  vector  (radians) 

c dc  is  the  star  declination  vector 

c rv  is  the  radial  velocity  vector 

c pmra  is  the  right  ascension  proper  motion  vector 
c pmdc  is  the  declination  proper  motion  vector 

c m3  is  the  number  of  condition  equations 

c 

C OUTPUTS : 

c f is  the  vector  of  the  condition  equations  evaluated  at  the 

current 

c values  of  the  observations  and  the  parameters 

c 

c declare  and  dimension  variables 

c 

implicit  real*8  (a-h,o-z) 

dimension  param (*) , ra (*) , dc ( *) , rv ( *) , pmra (*) ,pmdc (*) , f (m3) 
c 

c start  loop 

c 
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do  10  i=l, m3, 3 
knt= (i-1) /3+1 
i2=i+l 
i3=i+2 

f (i)  =rv(knt)  -param(l)  *dcos  (ra  (knt)  ) *dcos  (dc  (knt)  ) 

& -param(2) *dsin (ra (knt) ) *dcos (dc (knt) ) 

& -param(3) *dsin (dc (knt) ) 

f ( 12 ) =pmra (knt) +param(4) * (param(l) *dsin (ra (knt) ) 

& -param(2) *dcos (ra (knt) ) ) 

f ( 13 ) =pmdc (knt) +param(4) * (param(l) *dcos (ra (knt) ) *dsin (dc (knt) ) 
& +param(2) *dsin (ra (knt) ) *dsin (dc (knt) ) 

& -param(3) *dcos (dc (knt) ) ) 

10  continue 
c 

c finished 

c 

return 

end 

c 

subroutine  getsig (epmra, epmdc, erv, eplx, nobsv,m,ms, sigma) 
c 

c INPUTS : 

c epmra  is  the  vector  of  errors  in  the  ra  proper  motions 

c epmdc  is  the  vector  of  errors  in  the  dec  proper  motions 

c erv  is  the  vector  of  errors  in  the  radial  velocities 
c eplx  is  the  vector  of  errors  in  the  parallaxes 
c m is  the  number  of  observations 

c ms  is  the  number  of  rows  (or  columns)  of  sigma 
c 

c OUTPUTS : 

c sigma  is  the  covariance  matrix  of  the  observations 

c 

c declare  and  dimension  variables 

c 

implicit  real*8  (a-h,o-z) 

dimension  epmra (*) ,epmdc(*) ,erv(*) ,eplx(*) ,nobsv(*) 
dimension  sigma (ms, ms) 
c 

c start  loop 

c 

jpos=l 
do  10  i=l,m 
j2=jpos+l 
j3=jpos+2 
j4= jpos+3 

sigma ( jpos, jpos) =erv(i) *erv(i) 
sigma ( j 2 , j 2 ) =epmra (i) *epmra (i) 
sigma ( j3, j3) =epmdc (i) *epmdc (i) 
sigma ( j 4 , j 4) =eplx (i) *eplx (i) 
jpos= jpos+nobsv(i) 

10  continue 
return 
end 
c 

subroutine  getrho (sxdot, sydot, szdot, splx,m4, rho) 
c 

C INPUTS: 

c sxdot  is  the  velocity  constraint  on  the  x component 

c sydot  is  the  velocity  constraint  on  the  y component 

c szdot  is  the  velocity  constraint  on  the  z component 
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c splx  is  the  parallax  constraint,  obtained  from  cluster  angular 
c diameter 

c m4  is  the  row  or  column  dimension  of  rho 
c 

c OUTPUTS : 

c rho  is  the  covariance  matrix  of  the  constrained  parameters 

c 

c declare  and  dimension  variables 

c 

implicit  real*8  (a-h,o-z) 
dimension  rho(m4,m4) 
c 

c start  loop 

c 

do  10  i=l,m4,4 
i2=i+l 
i3=i+2 
i4=i+3 

rho (i,  i) =sxdot*sxdot 
rho (i2, i2) =sydot*sydot 
rho (i3, i3) =szdot*szdot 
rho (i4, i4) =splx*splx 
10  continue 
c 

c finished 

c 

return 

end 

c 

subroutine  getb (a, m3, n,m4, nobsv, b) 
c 

C INPUTS : 

c a is  the  jacobian  matrix  of  partials  with  respect  to  the 

parameters 

c m3,  n are  the  row  and  column  dimensions  of  a 

c m4  is  the  column  dimension  of  b 

c nobsv  is  the  vector  containing  the  number  of  observations  for  each 

c star  (2,  3 or  4) 

c 

C OUTPUTS : 

c b is  the  Jacobian  matrix  of  the  condition  equations  with  respect 

to  the 

c constrained  parameters 

c 

c declare  and  dimension  variables 

c 

implicit  real*8  (a-h,o-z) 
dimension  a (m3, n) , nobsv (*) ,b (m3,m4) 
c 

c start  loop 

c 

j=0 

do  20  i=l,m3,3 
i2=i+l 
i3=i+2 
do  10  k=l, 4 
ic= j+k 

b (i,  ic)  =a  (i,  k) 
b (i2,  ic) =a  (i2,k) 
b (i3,  ic)  =a  (i3,  k) 
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10  continue 
j=j+4 
20  continue 
c 

c finished 

c 

return 

end 

c 

subroutine  getdta (ra,dc, pmra, epmra, pmdc, epmdc, rv, erv,plx,eplx, 

& nobsv,m) 

c 

c INPUTS : 

c m is  the  number  of  stars 

c 

c OUTPUTS : 

c ra  is  the  vector  of  star  ra's,  units  are  radians 

c dc  is  the  vector  of  star  declinations,  units  are  radians 

c pmra  is  the  vector  of  star  ra  proper  motions,  units  are  radians 
c epmra  is  the  vector  of  standard  errors  associated  with  pmra 

c pmdc  is  the  vector  of  star  declination  proper  motions,  units  are 

radians 

c epmdc  is  the  vector  of  standard  errors  associated  with  pmdc 

c rv  is  the  vector  of  star  radial  velocities,  units  are  pc  per  yr 

c erv  is  the  vector  of  standard  errors  associated  with  rv 

c plx  is  the  vector  of  star  parallaxes,  units  are  seconds 

c eplx  is  the  vector  of  standard  errors  associated  with  plx 
c nobsv  is  the  vector  containing  the  number  of  observations  for  each 
c star  (2,  3 or  4) 

c 

c declare  and  dimension  variables 

c 

implicit  real*8  (a-h,o-z) 
real*8  ks2pyr 

dimension  ra  (*) ,dc  (*) ,pmra (*) ,pmdc (*) , rv(*) ,plx (*) 
dimension  epmra (*) , epmdc (*) , erv (*) , eplx (*) , nobsv (*) 
c 

c establish  constants  and  conversions 

c 

pi=3 . 14159654d0 
dg2rad=pi/l . 8d2 
ks2pyr=l . 027712166D-6 
sc2rad=dg2rad/3 . 6d3 
c 

c open  file 

c 

open  (4,  file=' stars. in') 
c 

c start  loop 

c 

do  10  i=l,m 

read  (4, 400, end=90)  alpha, delta, pra, epra, pdc, epdc, plx(i) , 

& eplx (i) , rho, erho, nobsv(i) 

400  format (3x,  2fll.6,6f9.5,2f6.2,i2) 

c 

c convert  to  the  appropriate  units 

c 

ra  (i) =alpha*dg2rad 
dc (i) =delta*dg2rad 
pmra (i) =pra*sc2rad 
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epmra (i) =epra*sc2rad 
pmdc (i) =pdc*sc2rad 
epmdc (i) =epdc*sc2rad 
rv (i) =rho*ks2pyr 
erv (i) =erho*ks2pyr 
continue 

finished 

close  (4) 

return 

end 

subroutine  spsinv (array, m3) 

INPUTS: 

array  is  a square  sparse  matrix  of  dimension  m3  with  3x3  blocks 
m3  is  dimension  of  array 

OUTPUTS: 

array  is  returned  as  its  inverse  (inverted  in  place) 

dimension  and  declare  variables 

implicit  real*8  (a-h,o-z) 
dimension  array (m3, m3) , temp(3,3) 

start  outer  loop 


jpos=0 

do  20  i=l,m3,3 
i2=i+l 
i3=i+2 


put  block  into  temporary  array 

do  5 j=l,3 
jl=j+jpos 

temp (1, j) =array (i, jl) 
temp (2, j)=array(i2, jl) 
temp ( 3 , j ) =array ( i3 , j 1 ) 

invert  temporary  array 

call  syminv(temp, 3) 

put  inverted  block  into  place 

do  10  j=l,3 
jl=j+jpos 

array (i, jl)=temp(l, j) 
array (i2, jl) =temp (2, j) 
array(i3, jl)=temp(3, j) 

increment  column 

jpos= jpos+3 


continue 
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c finished 
c 

return  * 

end 

subroutine  ef fcal (evec, cvec, n, ef fey) 
c 

c this  subroutine  will  calculate  the  'efficiency'  of  a fit,  given 
c the  vector  of  the  diagonal  of  the  scaled  covariance  matrix,  and 

c its  associated  eigenvalues 

c 

implicit  real*8  (a-h,o-z) 
dimension  evec (n) , cvec (n) 
c 

c calculate  the  product  of  the  eigenvalues 

c 

eprd=l . OdO 
do  10  i=l,n 
10  eprd=evec (i) *eprd 
c 

c calculate  the  product  of  the  covariance  diagonal  elements 
c 

cprd=l . OdO 
do  20  i=l,n 
20  cprd=cvec (i) *cprd 
c 

c the  efficiency  is  now  equal  to  eprd/cprd  to  the  1/nth  power 
c 

epwr=l . OdO/n 
ef fcy= (eprd/cprd) **epwr 
c 

c finished 

c 

return 

end 

C234567 

subroutine  clearv (ndim, vector) 
c 

c this  subroutine  will  fill  a vector  with  zeros 

c 

implicit  real*8  (a-h,o-z) 
dimension  vector (*) 
c 

c fill  vector  with  0's 

c 

do  10  i=l,ndim 
10  vector (i) =0 . OdO 
c 

c finished 

c 

return 

end 

c 

C234567 

subroutine  clearr (nrow, ncol, array) 
c 

c this  subroutine  will  fill  an  array  with  zeros 

c 

implicit  real*8  (a-h,o-z) 
dimension  array (nrow, ncol) 
c 


onoooo  ooo  ooo  oo 
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c fill  array  with  0's 

c 

do  10  i=l,nrow 
do  10  j=l,ncol 
10  array (i, j) =0 . OdO 
c 

c finished 

c 

return 

end 

c 

SUBROUTINE  VTXVG  (VI, V2 , NDIM, VALUE) 

This  subroutine  multiplies  a row  vector  by  a column  vector, 
outputting  a scalar,  VALUE 
IMPLICIT  REAL* 8 (A-H,0-Z) 

DIMENSION  VI (*) , V2 (*) 

VALUE=0 . ODO 
DO  10  1=1, NDIM 
1 0 VALUE=V1 ( I ) * V2 ( I ) +VALUE 
RETURN 
END 

SUBROUTINE  MATADD  ( ARRAY1 , ARRAY2 , NROW, NCOL, ARRAYO) 

This  subroutine  will  add  2 double  precision  arrays  of 
arbitrary  size 
IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  ARRAY 1 (NROW, NCOL) , ARRAY 2 (NROW, NCOL) 

DIMENSION  ARRAYO (NROW, NCOL) 

DO  10  1=1, NROW 
DO  10  J=l,NCOL 

1 0 ARRAYO ( I , J) =ARRAY1 ( I , J) + ARRAY 2 ( I , J) 

RETURN 
END 

SUBROUTINE  MATSUB  (ARRAY1, ARRAY2, NROW, NCOL, ARRAYO) 

This  subroutine  will  subtract  two  double  precision  arrays 
of  arbitrary  size 
IMPLICIT  REAL* 8 (A-H,0-Z) 

DIMENSION  ARRAY 1 (NROW, NCOL) , ARRAY2 (NROW, NCOL) 

DIMENSION  ARRAYO (NROW, NCOL) 

DO  10  1=1, NROW 
DO  10  J=l, NCOL 

10  ARRAYO (I, J)=ARRAY1 (I, J) -ARRAY 2 (I, J) 

RETURN 
END 

SUBROUTINE  NEGATE  (NROW, VECTOR) 

This  subroutine  returns  the  negative  of  VECTOR.  VECTOR  is  replaced 
by  its  negative. 

Dimension  array 
IMPLICIT  REAL* 8 (A-H,0-Z) 

DIMENSION  VECTOR (*) 

C Negate  VECTOR 

DO  10  1=1, NROW 
10  VECTOR (I) =-1.0* VECTOR (I) 

C Finished 

RETURN 
END 


non 


c 
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SUBROUTINE  MTXMG  (Ml,  M2,  NCI,  NR1R2 , NC2,  MOUT) 
DOUBLE  PRECISION  Ml (NR1R2,NC1) , M2 (NR1R2 , NC2 ) 
INTEGER* 4 NCI,  NR1R2 , NC2 
DOUBLE  PRECISION  MOUT (NCI, NC2) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c- 


Multiply  the  transpose  of  a matrix  and  a matrix,  of  arbitrary  size. 


VARIABLE  I/O 


Ml  I 

M2  I 

NCI  I 

NR1R2  I 

NC2  I 

MOUT  0 


DESCRIPTION 

Left-hand  matrix  whose  transpose  is  to  be  multiplied 
Right-hand  matrix  to  be  multiplied. 

Column  dimension  of  Ml  and  row  dimension  of  MOUT. 

Row  dimension  of  Ml  and  row  dimension  of  M2 . 

Column  dimension  of  M2  and  column  dimension  of  MOUT. 
Product  matrix  M1**T  * M2. 

MOUT  must  not  overwrite  either  Ml  or  M2 . 


VERSION  1.00  — 1984  Nov  08  — W.M.  Owen 


DOUBLE  PRECISION  SUM 
INTEGER  I, J, K 


C 

C Perform  the  matrix  multiplication 
C 

DO  10  1=1, NCI 
DO  10  J=1,NC2 
SUM  = 0 .DO 
DO  5 K=1,NR1R2 

5 SUM  = SUM  + Ml  (K,  I)  *M2  (K,  J) 

10  MOUT  (I,  J)  = SUM 

C 

RETURN 

END 


SUBROUTINE  MTXVG  (Ml,  V2,  NCI,  NR1R2 , VOUT) 
DOUBLE  PRECISION  Ml (NR1R2 , NCI ) , V2(NR1R2) 
INTEGER* 4 NCI,  NR1R2 
DOUBLE  PRECISION  VOUT (NCI) 


C 

C 

C 


Multiply  the  transpose  of  a matrix  and  a vector  of  arbitrary  size. 


c 

VARIABLE 

I/O 

DESCRIPTION 

c 

Ml 

I 

Left-hand  matrix  whose  transpose  is  to  be  multiplied 

c 

V2 

I 

Right-hand  vector  to  be  multiplied. 

c 

NCI 

I 

Column  dimension  of  Ml  and  length  of 

VOUT. 

c 

NR1R2 

I 

Row  dimension  of  Ml  and  length  of  V2 

• 

c 

c 

VOUT 

O 

Product  vector  M1**T  * V2 . 

VOUT  must  not  overwrite  either  Ml  or 

V2  . 

C 

C 

C 

C- 


VERSION  1.00  — 1984  Nov  08  — W.M.  Owen 


DOUBLE  PRECISION  SUM 
INTEGER  I, K 


Perform  the  matrix-vector  multiplication 


DO  10  1=1, NCI 
SUM  = 0 .DO 


non 
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c 


c 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 


DO  5 K=1 , NR1R2 

5 SUM  = SUM  + Ml (K, I) *V2 (K) 

10  VOUT(I)  = SUM 

RETURN 

END 


SUBROUTINE  MXMG  (Ml,  M2,  NR1,  NC1R2,  NC2,  MOUT) 
DOUBLE  PRECISION  Ml (NR1 , NC1R2 ) , M2 (NC1R2,NC2) 
INTEGER* 4 NR1,  NC1R2,  NC2 
DOUBLE  PRECISION  MOUT (NR1,NC2) 


Multiply  two  matrices  of  arbitrary  size. 


VARIABLE  I/O 


Ml  I 

M2  I 

NR1  I 

NC1R2  I 

NC2  I 

MOUT  0 


DESCRIPTION 

Left-hand  matrix  to  be  multiplied. 

Right-hand  matrix  to  be  multiplied. 

Row  dimension  of  Ml  (and  also  MOUT) . 

Column  dimension  of  Ml  and  row  dimension  of  M2 . 
Column  dimension  of  M2  (and  also  MOUT) . 

Product  matrix  M1*M2. 

MOUT  must  not  overwrite  either  Ml  or  M2 . 


VERSION  1.00  — 1984  Nov  08  — W.M.  Owen 


DOUBLE  PRECISION  SUM 
INTEGER  I , J, K 


Perform  the  matrix  multiplication 

DO  10  1=1, NR1 
DO  10  J=1 , NC2 
SUM  = 0 .DO 
DO  5 K=l, NC1R2 

5 SUM  = SUM  + Ml (I,K) *M2 (K, J) 

10  MOUT (I, J)  = SUM 

RETURN 

END 


SUBROUTINE  MXMTG  (Ml,  M2,  NR1,  NC1C2 , NR2 , MOUT) 

DOUBLE  PRECISION  Ml (NR1,NC1C2) , M2 (NR2 , NC1C2 ) 

INTEGER* 4 NR1,  NC1C2 , NR2 
DOUBLE  PRECISION  MOUT (NR1, NR2 ) 

Multiply  a matrix  and  the  transpose  of  a matrix,  of  arbitrary  size 


VARIABLE  I/O 


Ml  I 

M2  I 

NR1  I 

NC1C2  I 

NR2  I 

MOUT  O 


DESCRIPTION 

Left-hand  matrix  to  be  multiplied. 

Right-hand  matrix  whose  transpose  is  to  be  mult. 
Row  dimension  of  Ml  and  row  dimension  of  MOUT. 
Column  dimension  of  Ml  and  column  dimension  of  M2 
Row  dimension  of  M2  and  column  dimension  of  MOUT. 
Product  matrix  Ml  * M2**T. 

MOUT  must  not  overwrite  either  Ml  or  M2 . 


VERSION  1.01  — 1984  Nov  12  — W.M.  Owen 


ooo  ooonooonooooon  n o ooo 
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DOUBLE  PRECISION  SUM 
INTEGER  I,J,K 

Perform  the  matrix  multiplication 

DO  10  1=1, NR1 
DO  10  J=l, NR2 
SUM  = 0 .DO 
DO  5 K=l, NC1C2 

5 SUM  = SUM  + Ml  (I,  K)  *M2  ( J,  K) 

10  MOUT (I, J)  = SUM 


RETURN 

END 

SUBROUTINE  MXVG  (Ml,  V2,  NR1,  NC1R2 , VOUT) 
DOUBLE  PRECISION  Ml (NR1,NC1R2) , V2(NC1R2) 
INTEGER* 4 NR1,  NC1R2 
DOUBLE  PRECISION  VOUT(NRl) 

Multiply  a matrix  and  a vector  of  arbitrary  size. 


VARIABLE  I/O 


Ml  I 

V2  I 

NR1  I 

NC1R2  I 

VOUT  O 


DESCRIPTION 

Left-hand  matrix  to  be  multiplied. 
Right-hand  vector  to  be  multiplied. 

Row  dimension  of  Ml  and  length  of  VOUT. 
Column  dimension  of  Ml  and  length  of  V2 
Product  vector  M1*V2. 

VOUT  must  not  overwrite  either  Ml  or  V2 


VERSION  1.00  — 1984  Nov  08  — W.M.  Owen 


DOUBLE  PRECISION  SUM 
INTEGER  I, K 

Perform  the  matrix-vector  multiplication 

DO  10  1=1, NR1 
SUM  = 0 .DO 
DO  5 K=l, NC1R2 

5 SUM  = SUM  + Ml (I,K) *V2 (K) 

10  VOUT (I)  = SUM 
C 

RETURN 

END 

c 

SUBROUTINE  VADDG  (VI,  V2,  NDIM,  VOUT) 
DOUBLE  PRECISION  VI (NDIM) , V2 (NDIM) 
INTEGER* 4 NDIM 
DOUBLE  PRECISION  VOUT (NDIM) 

C 


c 

n 

Add  two 

vectors 

of  arbitrary  size. 

c 

VARIABLE 

I/O 

DESCRIPTION 

c 

VI 

I 

First  vector  to  be  added. 

c 

V2 

I 

Second  vector  to  be  added. 

c 

NDIM 

I 

Dimension  of  VI,  V2,  and  VOUT. 

c 

c 

VOUT 

O 

Sum  vector,  VI  + V2 . 

VOUT  can  overwrite  either  VI  or  V2 

o o o n 
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VERSION  1.00  — 1984  Nov  08  — W.M.  Owen 


INTEGER  I 
C 

DO  1 1=1 , NDIM 

1 VOUT(I)  = VI (I)  + V2 ( I ) 

C 

RETURN 
END 
c 

SUBROUTINE  VSUBG  (VI,  V2,  NDIM,  VOUT) 

DOUBLE  PRECISION  VI (NDIM) , V2 (NDIM) 

INTEGER* 4 NDIM 
DOUBLE  PRECISION  VOUT (NDIM) 

C 

C Subtract  two  vectors  of  arbitrary  size. 

C 

C VARIABLE 
C VI 
C V2 
C NDIM 
C VOUT 
C 
C 

C VERSION  1.00  — 1984  Nov  08  — W.M.  Owen 
C 

C- 

INTEGER  I 
C 

DO  1 1=1, NDIM 

1 VOUT  (I)  = VI  (I)  - V2  (I) 

C 

RETURN 

END 

c 

SUBROUTINE  SYMINV (ARRAY, NDIM) 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  ARRAY (NDIM, NDIM) , S (50, 50) , ARYINV(50, 50) 
C CALCULATE  S ( 1, 1) , S (1, J) 

S(1,1)=DSQRT (ARRAY (1,1) ) 

DO  10  J=2,NDIM 
10  S(1,J)=ARRAY(1,J)/S(1,1) 

C CALCULATE  S (I, I) , S (I, J) 

1=2 

SUM3=0 . 0D0 
15  SUM2=0 . 0D0 
LOOPLM=I-l 
DO  20  L=1 , LOOPLM 
20  SUM2=S (L, I) *S (L, I) +SUM2 

S (I,  I) =DSQRT (ARRAY (I, I) -SUM2) 

J=I+1 

25  DO  30  L=l, LOOPLM 
30  SUM3=S (L, I) *S (L,J) +SUM3 

S (I,  J)  = (ARRAY  (I,  J)  -SUM3)  /S  (I,  I) 

J=J+1 

SUM3=0 . 0D0 

IF  (J  .LE.  NDIM)  GOTO  25 
1=1+1 


I/O  DESCRIPTION 

I First  vector  (minuend) . 

I Second  vector  (subtrahend) . 

I Dimension  of  VI,  V2,  and  VOUT. 

O Difference  vector,  VI  - V2 . 

VOUT  can  overwrite  either  VI  or  V2 . 
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IF  (I  .LE.  NDIM)  GOTO  15 
C CALCULATE  INVERSE 

CD=0 . 0D0 
SUM=0 . 0D0 

ARYINV  (NDIM,  NDIM)  =1 .0D0/  ( S (NDIM,  NDIM)  *S  (NDIM,  NDIM)  ) 

I=NDIM-1 
40  K=ND IM 

45  IF  (K  .EQ.  I)  CD=1 . ODO/S (1,1) 

IF  (K  .NE.  I)  CD=0 . 0D0 
LOOPST=I+l 

DO  50  NU=LOOP ST , ND IM 
50  SUM=S ( I , NU) *ARYINV (NU, K) +SUM 
ARYINV ( I , K)  = (CD-SUM) /S (I,  I) 

ARYINV (K, I ) = ARYINV ( I , K) 

SUM=0 . 0D0 
K=K-1 

IF  (K  .GE.  I)  GOTO  45 
1=1-1 

IF  (I  .GE.  1)  GOTO  40 
DO  55  12=1, NDIM 
DO  55  J2=l, NDIM 
55  ARRAY (12, J2) =ARYINV(I2, J2) 

RETURN 

END 

c 

C234567 

subroutine  jacobi (a,n,np,d,v,nrot) 
c 

c this  subroutine  computes  the  eigenvectors  and  eigenvalues  of  a 
real 

c symmetric  matrix  a,  of  size  n by  n,  stored  in  a physical  np  by  np 

c array.  On  output  the  elements  of  a above  the  diagonal  are 

destroyed. 

c d returns  the  eigenvalues  of  a in  its  first  n elements,  v is  a 
matrix 

c with  the  same  logical  and  physical  dimensions  as  a whose  columns 

c contain,  on  output,  the  normalized  eigenvectors  of  a.  nrot  returns 

c the  number  of  required  jacobi  rotations, 
c 

c this  program  assumes  underflows  are  treated  as  zeros!!!!! 
c 

implicit  real*8  (a-h,o-z) 
parameter  (nmax=100) 

dimension  a (np,np) ,d(np) , v(np,np) ,b (nmax) , z (nmax) 
do  12  ip=l,n 
do  11  iq=l,n 

v (ip, iq) =0 . OdO 

11  continue 
v(ip, ip) =1 . OdO 

12  continue 
do  13  ip=l,n 

b (ip)  =a  (ip,  ip) 
d (ip)  =b  (ip) 
z (ip) =0 . OdO 
continue 
nrot=0 

do  24  i=l, 50 
sm=0 . OdO 
do  15  ip=l, n-1 

do  14  iq=ip+l,n 


13 
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sm=sm+dabs (a (ip, iq) ) 
continue 
continue 

if  (sm.eq. 0 . OdO)  return 
if  (i.lt.4)  then 

tresh=0 .2*sm/n**2 
else 

tresh=0 . OdO 
endif 

do  22  ip=l,n-l 

do  21  iq=ip+l,n 

g=l . 0d2*dabs (a (ip, iq) ) 

if  ( (i.gt .4) .and. (dabs (d(ip) ) tg.eq.dabs (d(ip) ) ) 
.and. (dabs (d(iq) ) tg.eq.dabs (d(iq) ) ) ) then 
a (ip,  iq) =0 . OdO 

else  if  (dabs (a (ip, iq) ) .gt .tresh)  then 
h=d (iq) -d(ip) 

if  (dabs (h) +g.eq. dabs (h) ) then 
t=a (ip, iq) /h 
else 

theta=0 . 5*h/a (ip, iq) 

t=l . OdO/ (dabs (theta) +dsqrt (1 . 0d0+theta**2) ) 
if  (theta. It. 0. OdO)  t=-t 
endif 

c=l . OdO/dsqrt (l+t**2) 
s=t*c 

tau=s/ (l.OdO+c) 
h=t*a (ip, iq) 
z (ip)  =z  (ip)  -h 
z (iq)  =z  (iq)  +h 
d(ip)  =d(ip)  -h 
d(iq) =d (iq) +h 
a (ip, iq) =0 . OdO 
do  16  j=l,ip-l 
g=a ( j , ip) 
h=a ( j, iq) 

a ( j, ip) =g-s* (h+g*tau) 
a ( j , iq) =h+s* (g-h*tau) 
continue 

do  17  j=ip+l, iq-1 
g=a (ip, j) 
h=a ( j, iq) 

a (ip, j) =g-s* (h+g*tau) 
a ( j, iq) =h+s* (g-h*tau) 
continue 
do  18  j=iq+l,n 
g=a (ip, j ) 
h=a (iq, j) 

a (ip, j ) =g-s* (h+g*tau) 
a (iq, j)=h+s* (g-h*tau) 
continue 
do  19  j=l,n 
g=v( j,ip) 
h=v( j, iq) 

v( j,ip)=g-s* (h+g*tau) 
v( j, iq) =h+s* (g-h*tau) 
continue 
nrot=nrot+l 
endif 
continue 
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22  continue 

do  23  ip=l,n 

b (ip)  =b  (ip)  +z  (ip) 
d(ip) =b (ip) 
z (ip) =0 . OdO 

23  continue 

24  continue 

write (*,  *)  '50  iterations  should  never  happen!' 

return 
end 
c 

subroutine  eigsrt (d, v, n, np) 
c 

c this  subroutine  sorts  the  eigenvalues,  d,  and  the  eigenvectors,  v, 

c as  output  by  the  jacobi  routine.  It  will  put  the  eigenvalues  into 

c ascending  order,  and  then  rearrange  the  eigenvectors 

correspondingly . 
c 

implicit  real*8  (a-h,o-z) 
dimension  d (np) , v (np, np) 
do  13  i=l, n-1 
k=i 
p=d(i) 

do  11  j=i+l,n 

if  (d(j) .ge.p)  then 
k=j 
p=d ( j ) 
endif 

11  continue 

if  (k.ne.i)  then 
d (k) =d(i) 
d(i)=p 
do  12  j=l,n 
p=v( j, i) 
v(  j,  i)  =v(  j,k) 
v(  j,k)=p 

12  continue 
endif 

13  continue 
return 
end 

C234567 

subroutine  datard  (bdnum,  ra,dc,pmra,epmra,pmdc,epmdc,  rv, 

& erv, plx, eplx, nobsv, vmag, blessv, m) 

c 

c this  subroutine  reads  the  data  file  schwan.dat,  converting 
c everything  to  the  appropriate  units... 

c 

implicit  real*8  (a-h,o-z) 

dimension  bdnum (*) ,ra (*) ,dc(*) ,pmra (*) ,pmdc(*) ,rv(*) 
dimension  epmra ( *) , epmdc ( *) , erv ( *) , vmag ( *) , blessv ( *) 
dimension  nobsv (*) ,plx(*) ,eplx(*) 
character*l  indrv 
character*7  hdnum 
c 

c establish  pi,  and  seconds  to  radians,  km/s  to  pc/yr  conversions 
c 

pi=3.141592654d0 
sc2rad=pi/ (1 . 8d2*3 . 6d3) 
rks2py=l . 022712166d-6 


Ill 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


open  data  file 

open  (1, f ile= ' schwan2 .dat ' ) 

read  data 


knt=l 

5 read(l, 100, end=99) 
& 

& 


ivbnum, bdnum(knt) ,hdnum, ra (knt) ,dc (knt) , 
ipmra, iepmra, ipmdc, iepmdc, rv (knt) , indrv, 
vmag (knt) , blessv (knt) , plx (knt) , ieplx 


convert  everything  to  the  appropriate  units... 
first  the  proper  motion  errors  to  radians/year 


if  (iepmra .eq. 0)  iepmra=l 
if  ( iepmdc . eq . 0 ) iepmdc=l 
epmra (knt) =iepmra*sc2rad/l . 0d4 
epmdc (knt) =iepmdc*sc2rad/l . 0d4 

get  radial  velocity  error,  convert  to  pc/yr 
if  (indrv. eq. 'G' ) erv (knt ) =0 . 5*rks2py 
if  (indrv. eq. 'W' ) erv(knt) =5 . 0*rks2py 

correct  the  radial  velocities  for  zero  point  error... 


rv (knt) =rv (knt) *rks2py 

convert  the  proper  motions  to  radians  per  year. . . 

pmra (knt) =ipmra/l . 0d4*sc2rad 
pmdc (knt) =ipmdc/l . 0d4*sc2rad 


get  the  parallax  errors 
eplx (knt) =ieplx/l . 0d4 

store  number  of  observations  for  this  star 


nobsv(knt) =3 

if  (ieplx  .gt.  0)  nobsv(knt)=4 

everything  is  now  in  the  right  units,  so  increment 
the  counter  and  get  another  star. 

knt=knt+l 
goto  5 

format  statements 

100  format  (i3,  f 8 .3,  a7, 2fl2 .8,  lx,  i4,  i3,  lx,  i4,  i3,  f5 . 1,  lx,  al,  f 6 .2,  f5 .2, 
& f 8 . 4, i4) 


finished,  establish  number  of  stars. 


99  m=knt-l 
close  (1) 
return 
end 


c 

c 

c 


APPENDIX  B 

TABLE  OF  ADJUSTMENT  RESIDUALS 


Table  B-l.  Observation  residuals. 


No. 

v„r  (km  s1) 

v„„  ("  yr1) 

V„s  ("  yr1) 

v,  (arcsec.) 

l 

0.1110 

0.0000 

0.0000 

0.0080 

2 

-0.0180 

0.0000 

-2.0000e-5 

0.0002 

3 

0.3230 

0.0000 

1.0000e-5 

0.0096 

4 

0.5660 

1.0000e-5 

-2.0000e-5 

-0.0101 

5 

-0.7360 

-0.0001 

-2.0000e-5 

6 

1.0890 

3.0000e-5 

0.0003 

0.0184 

7 

3.5680 

4.0000e-5 

0.0003 

0.0045 

8 

3.7390 

0.0000 

-4.0000e-5 

-0.0097 

9 

-1.6300 

0.0004 

0.0002 

10 

7.3640 

1.0000e-5 

-0.0002 

-0.0258 

11 

2.2500 

-0.0001 

0.0006 

0.0377 

12 

3.1250 

0.0001 

0.0001 

13 

4.0410 

-0.0002 

-0.0005 

0.0100 

14 

7.2350 

-6.0000e-5 

-0.0006 

15 

-3.1980 

-0.0003 

-0.0018 

0.0063 

16 

5.6460 

0.0002 

0.0006 

17 

4.5900 

-0.0001 

-0.0015 

18 

1.7510 

0.0001 

-0.0029 

-0.0409 

19 

-2.1630 

-0.0008 

-0.0001 

20 

2.1530 

0.0002 

0.0011 

21 

-2.2790 

-0.0004 

-0.0009 

22 

3.4160 

-9.0000e-5 

-0.0008 

23 

2.7620 

0.0003 

0.0011 

24 

-1.2620 

0.0002 

0.0020 

25 

6.6300 

-0.0050 

-0.0045 

26 

4.7360 

-0.0008 

-0.0049 

27 

2.2690 

0.0001 

0.0006 

0.0044 

28 

3.0370 

-6.0000e-5 

-7.0000e-5 

29 

-0.8880 

-0.0007 

-0.0016 

0.0030 

30 

7.1150 

-0.0007 

-0.0011 

31 

-0.1520 

0.0001 

0.0014 

0.0022 

32 

3.6840 

0.0002 

0.0009 

-0.0078 

33 

3.3430 

4.0000e-5 

-1.0000e-5 

34 

0.1270 

0.0002 

0.0011 

0.0045 

35 

0.9590 

-6.0000e-5 

-0.0003 

0.0068 

36 

3.0950 

-0.0003 

-0.0014 

0.0035 

37 

1.3840 

-3.0000e-5 

-0.0002 

38 

2.0520 

-0.0003 

-0.0005 

39 

4.7420 

0.0003 

0.0017 
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Table  B-l  — continued 


No. 

v„r  (km  s-1) 

v„„  ("  yr1) 

v..s(  y^1) 

40 

-0.2640 

2.0000e-5 

0.0012 

41 

2.9890 

0.0004 

0.0014 

42 

0.2380 

0.0004 

0.0015 

43 

3.0100 

-0.0002 

-0.0005 

44 

-0.6800 

4.0000e-5 

0.0006 

45 

-0.0640 

-0.0032 

-0.0050 

46 

-0.0250 

-0.0012 

6.0000e-5 

47 

-0.4030 

-0.0012 

-0.0044 

48 

-0.2030 

-0.0013 

0.0003 

49 

-0.0500 

-0.0006 

-0.0016 

50 

-0.0920 

0.0002 

0.0006 

51 

1.8790 

-0.0004 

-0.0011 

52 

-0.0250 

0.0006 

0.0007 

53 

4.4300 

0.0004 

2.0000e-5 

54 

0.1390 

-0.0004 

0.0003 

55 

0.3260 

-0.0013 

-0.0031 

56 

-0.1480 

3.0000e-5 

-0.0007 

57 

-0.2090 

-0.0010 

-0.0028 

58 

-3.1370 

-0.0016 

-0.0021 

59 

-0.4130 

0.0021 

-0.0034 

60 

-0.0620 

2.0000e-5 

-2.0000e-5 

61 

0.1760 

-0.0024 

-0.0019 

62 

-1.6300 

0.0030 

-0.0039 

63 

-0.2020 

-0.0004 

0.0004 

64 

-0.0380 

-0.0003 

-0.0004 

65 

0.8640 

-0.0027 

0.0017 

66 

5.0850 

-0.0010 

-0.0028 

67 

-0.1000 

5.0000e-5 

0.0024 

68 

-0.0720 

-0.0005 

-0.0041 

69 

0.3710 

0.0002 

0.0037 

70 

0.4670 

0.0001 

0.0001 

71 

-4.8990 

-0.0002 

-0.0007 

72 

-0.1430 

-0.0002 

0.0006 

73 

-0.3450 

0.0002 

0.0020 

74 

0.2200 

-0.00010 

-0.0041 

75 

0.1220 

-0.0003 

0.0013 

76 

0.1680 

-0.0004 

-0.0004 

77 

-2.7270 

-0.0009 

-0.0021 

78 

2.3620 

-0.0014 

-0.0043 

79 

0.3160 

-0.0004 

-0.0024 

80 

-0.3270 

-0.0002 

0.0019 

81 

-1.3810 

-0.0038 

-0.0052 

82 

-0.2360 

0.0022 

0.0050 

83 

0.3050 

0.0001 

0.0075 

84 

-0.1320 

-0.0008 

-0.0020 

85 

1.2870 

-0.0005 

-0.0011 

86 

0.0690 

0.0009 

-0.0001 

87 

-3.2740 

0.0011 

0.0011 

88 

0.5500 

-0.0028 

-0.0041 

89 

-0.4420 

-0.0019 

-0.0034 

(arc sec.) 


0.0064 


0.0236 

0.0037 


Table  B-l  — continued 
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No. 

v„r  (km  s1) 

V„„  ("  yr>) 

v,.s  ("  yr1) 

v,  (arc  sec.) 

90 

0.1160 

0.0004 

0.0021 

91 

-1.1300 

-0.0013 

-0.0017 

0.0054 

92 

0.3590 

-0.0009 

-0.0027 

0.0079 

93 

-0.0170 

-0.0003 

-0.0003 

0.0070 

94 

-0.2240 

0.0002 

0.0008 

95 

0.5400 

7.0000e-5 

0.0005 

0.0095 

96 

0.2520 

-0.0001 

-0.0011 

97 

-0.4950 

0.0019 

0.0019 

98 

-0.1490 

0.0005 

-0.0044 

99 

-3.5890 

-0.0004 

-0.0015 

100 

0.0330 

-2.0000e-5 

-0.0010 

-0.0012 

101 

0.1620 

-0.0001 

-0.0003 

0.0166 

102 

-4.4540 

0.0002 

-0.0018 

103 

0.2240 

0.0004 

0.0020 

104 

0.8560 

0.0004 

0.0018 

0.0073 

105 

-0.1850 

-0.0003 

-0.0002 

106 

-0.2370 

-3.0000e-5 

-0.0001 

107 

2.0950 

-0.0003 

-0.0014 

108 

-0.0530 

0.0005 

0.0002 

-0.0031 

109 

-2.8890 

0.0001 

0.0006 

110 

-0.1480 

0.0008 

0.0014 

-0.0031 

111 

-2.0710 

0.0026 

0.0082 

-0.0131 

112 

-0.3300 

0.0008 

0.0015 

113 

-0.2020 

0.0034 

0.0041 

-0.0144 

114 

6.1590 

-0.0002 

-0.0026 

-0.0057 

115 

3.4520 

0.0008 

0.0025 

116 

-0.0030 

0.0015 

0.0014 

117 

0.1640 

0.0011 

0.0036 

118 

-0.1160 

0.0039 

0.0078 

-0.0060 

119 

7.4920 

0.0018 

0.0038 

120 

-0.0940 

0.0003 

0.0012 

121 

-0.2930 

-0.0003 

-0.0010 

122 

0.3470 

0.0026 

0.0050 

123 

-0.0490 

0.0012 

0.0036 

124 

0.0470 

0.0010 

0.0028 

125 

0.2570 

0.0006 

0.0019 

126 

0.2630 

0.0004 

0.0007 

127 

0.7330 

-0.0004 

0.0007 

128 

0.2050 

0.0002 

-0.0005 

129 

-0.0460 

0.0018 

0.0036 

130 

0.3440 

0.0018 

0.0013 

131 

-0.0740 

0.0017 

0.0004 

132 

4.1070 

0.0024 

0.0011 

133 

7.3900 

0.0026 

-0.0005 

134 

-0.1200 

0.0011 

0.0002 

-0.0082 

135 

0.2090 

0.0031 

-0.0042 

136 

-0.5410 

-0.0005 

0.0054 

0.0068 

137 

-0.8690 

0.0008 

0.0019 

138 

0.9630 

-0.0003 

-0.0027 

139 

0.1310 

-0.0010 

-0.0012 

140 

0.7690 

0.0028 

0.0040 
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